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Chapter  1 
INTRODUCTION 


1 


l.l.  Statement  of  the  Problem 

The  purpose  of  this  study  is  to  incorporete  the  effects  of  multiple  boundary 
returns  end  non- isovelocity  conditions  into  a  time- varying  representation  of 
underwater  acoustic  reverberation  measured  from  a  moving  source.  The 
reverberation  model  used  as  a  foundation  for  this  study  is  outlined  by  Hodgkiss 
[1|.  The  model  was  origlnaily  developed  by  C.  L.  Ackerman  and  R.  L.  Kesser  at 
the  Applied  Research  Laboratory  of  The  Pennsylvania  State  University. 

Reverberation  can  be  described  as  a  lingering  of  sound  energy  in  a  space 
after  some  sound  has  been  generated.  In  the  ocean,  once  a  sound  has  been 
transmitted,  there  are  various  ways  in  which  the  energy  may  return  to  a  receiver 
situated  at  the  source  location.  The  ocean  environment  scatters  sound  in  all 
directions,  including  back  toward  the  direction  from  which  the  sound  originated. 
The  scattering  is  due  to  discontinuities  and  inhomogeneities  within  the  ocean. 
Scattering  functions  may  be  used  to  describe  the  backscattered  energy  levels 
as  a  function  of  time  and  frequency  for  different  portions  of  the  ocean.  The 
scattering  fonctioa  is  dependent  on  the  medium’s  scattering  strength,  i.e.,  its 
ability  to  scatter  sound. 

If  acoustic  beam  patterns  are  introduced  for  source  and  receiver  transducers, 
and  acoustic  spreading  and  absorption  are  considered,  a  scattering  function  of 
the  ocean  may  be  determined.  The  scattering  function  represents  a  complex, 
time-varying  filter,  through  which  a  transmitted  waveform  may  be  passed  to 
yield  the  reverberation  power  level  and  spectrum  versus  time  at  a  receiver. 

The  original  model  for  predicting  reverberation,  which  has  been  modified 
in  this  effort,  has  two  limitations.  First,  no  boundary  refiections  are  allowed. 
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The  ocean  surface  and  bottom  both  reflect  sound  forward  and  scalier  it  in 
all  directions  in  addition  to  scattering  the  sound  backward  toward  the  source. 
Secondly,  the  speed  of  sound  is  assumed  to  be  constant  throughout  the  ocean. 
In  the  real  ocean,  the  speed  of  sound  varies  from  point  to  point,  especially  over 
longer  distances. 

1.2.  Importance  of  the  Problem 

The  air- water  interface,  which  is  referred  to  as  the  surface  of  the  ocean, 
is  an  excellent  reflector  of  sound.  A  pressure-release  type  boundary  condition 
exists  and,  due  to  the  high  impedance  mismatch  between  the  water  and  the 
air,  the  impinging  energy  has  a  difficult  time  penetrating  the  boundary.  The 
consequence  is  that  most  of  the  energy  remains  contained  in  the  ocean,  either 
through  scattering  or  reflection. 

The  ocean  bottom  is  a  more  complex  boundary  than  the  surface.  In  some 
instances,  the  bottom  may  be  made  of  rock.  In  other  cases,  it  may  vary  with 
depth  from  soft  mud  to  harder  clay.  In  any  event,  reflection  does  occur  when 
sound  impinges  upon  the  ocean  bottom.  The  actual  amount  of  reflection  depends 
highly  on  the  composition  of  the  bottom.  Reflections  from  the  surface  and  the 
bottom  can  play  an  important  role  in  the  model,  as  sound,  which  is  backscattered 
after  the  occurrence  of  reflections,  may  contribute  significantly  to  the  received 
reverberation. 

In  the  actual  ocean,  the  speed  of  sound  is  by  no  means  constant.  Over 
somewhat  short  distances,  a  constant  speed  of  sound  may  be  a  good  assumption. 
However,  the  variation  in  the  sound  velocity  u  dependent  upon  such  factors  as 
temperature,  pressure,  and  salinity.  These  parameters  do  change  with  distance, 
mostly  with  depth.  Therefore,  a  non-isovelocity  medium  would  be  a  better 
assumption  in  a  realistic  model  of  the  ocean. 
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There  are  many  physical  phenomena  which  become  prevalent  in  a  non¬ 
isovelocity  environment.  With  a  varying  speed  of  sound,  spherical  spreading 
may  no  longer  be  assumed.  Acoustic  waves  au^e  refracted  and  their  paths  are 
no  longer  straight  lines.  Channels,  regions  in  which  a  portion  of  the  acoustic 
energy  becomes  trapped  between  two  horizontal  planes,  may  exist.  The  decay 
of  the  sound  pressure  level  over  distance  becomes  much  less  in  a  channel  than 
in  a  spherical-spreading  situation.  Convergence  zones  also  may  be  formed  in  a 
non-isovelocity  environment.  These  are  regions  in  which  the  sound  pressure  level 
is  greater  thsm  that  predicted  by  spherical  spreading.  In  contrast,  shadow  zone 
regions  in  which  the  sound  pressure  level  drops  drastically  over  distance,  could 
be  formed  in  an  environment  with  a  varying  sound  velocity.  All  of  these  factors 
add  to  the  complexity  of  predicting  underwater  reverberation  levels  and  spectral 
content. 

1.3.  Generalized  Approach 

The  reverberation  model  to  be  modified,  REVMOD  (REVerberation 
MODel),  consists  of  a  series  of  five  FORTRAN  computer  programs.  In  the  third 
program  of  the  series,  SCATTER,  a  scattering  function  at  a  particular  distance  is 
determined.  The  effects  that  the  sea  has  on  a  transmitted  signal  are  described  by 
the  scattering  function  [1].  The  carrier  frequency  and  pulse  length  are  necessary 
parameters,  but  the  signal  envelope  and  source  level  are  not  required  in  this 
section  of  REVMOD.  This  study  focuses  on  the  modification  of  SCATTER. 

In  SCATTER,  the  ocean  is  divided  into  a  spherical  grid  of  cells  at  a  given 
distance  from  the  source.  This  grid  has  a  thickness  proportional  to  the  pulse 
length.  Each  cell  contributes  only  one  type  of  reverberation,  either  surface, 
volume,  or  bottom  as  in  Figure  1.1.  In  Figure  1.1,  six  rows  of  ceils  are  shown. 
The  top  and  bottom  rows  consist  of  surface  and  bottom  cells,  respectively.  These 
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cells  are  two-dimensional,  and  are  therefore  represented  as  areas.  The  middle  four 
rows  are  volume  cells,  which  are  three  dimensional.  The  type  of  reverberation 
depends  on  where  the  cell  lies  in  the  ocean,  whether  it  is  on  the  surface,  in  the 
volume,  or  on  the  bottom.  The  grid  cells  are  processed  individually  for  their 
contributions  to  the  total  power  gain  of  the  reverberation.  The  term  “power 
gain”  is  used  because  the  actual  reverberation  power  level  is  not  determined 
until  a  transmitted  energy  level  is  specified.  This  energy  level  is  not  necessary  in 
the  determination  of  the  reverberation  power  gain.  For  each  cell,  a  backscattering 
cross-section  is  determined,  a  frequency  shift  due  to  the  velocities  of  both  the 
source  and  the  scatterer  is  introduced,  and  a  spreading  of  frequency  due  to  an 
assumed  Gaussian  distribution  of  random  velocity  at  the  scatterer  is  calculated. 
Acoustical  spreading  and  absorption  are  considered,  as  well  as  the  directivities 
of  the  source  and  the  receiver.  The  contributions  from  all  cells  are  summed  to 
obtain  the  power  gain  and  scattering  function  spectrum  at  a  given  distance. 

When  a  non-isovelocity  environment  and  multiple  forward  reflections  are 
considered,  straight-line  propagation  and  spherical  spreading  of  the  acoustic 
energy  no  longer  results.  In  order  to  solve  the  problem  of  non-spherical  spreading, 
the  concept  of  surfaces  of  constant  time  is  introduced.  With  straight-line 
propagation  and  no  boundary  reflections,  the  surface  of  constant  time  is  a  sphere. 
If  a  grid  similar  to  the  spherical  grid  in  SCATTER  is  to  be  utilized,  a  surface 
of  constant  time,  over  which  a  new  grid  is  set  up,  must  be  deflned.  Figure  1.2 
illustrates  a  series  of  constant  time  contours,  having  a  total  vertical  angular  span 
of  sixty  degrees,  and  centered  at  zero  degrees.  The  source  is  at  a  depth  of  flfty 
feet.  The  second  contour  in  this  series  shows  a  surface  reflection.  The  eighth 
shows  a  bottom  reflection  and  the  ninth  shows  a  surface-bottom  reflection.  Once 
this  new  grid  has  been  established,  other  ways  of  determining  the  backscattering 
cross  section,  Doppler  shift  and  frequency  spreading  must  be  formulated.  See 


<on«l«nl  tim«  contour* 


Comparison  of  a  spherically  spreading  wavefront 
downwardly  rdiacting  waveiront. 


8 


Figure  1.3  for  «  comp&rison  of  »  Msetioo  of  a  spoencaily  spreading  wavefront  for 
isovclocity  propagation,  the  sohd  curve,  to  the  same  section  of  a  downwardly 
refracting  wavefront,  the  dashed  curve,  which  is  due  to  a  decreasing  sound 
velocity  with  depth. 

Since  SCATTER  usee  a  constant  speed  of  sound  and  straight  line  propa¬ 
gation.  the  angle  at  which  a  ray  impinges  upon  a  scatterer  is  the  same  as  the 
transmitted  angle.  This  is  the  angle  used  for  Doppler  shift  calculations.  In  the 
non-isovelocity  case,  the  angle  at  the  scatterer  is,  in  general,  different  from  the 
transmitted  angle.  Knowledge  of  this  angle  is  very  important  when  the  scatterer 
is  moving  with  some  velocity.  This  velocity  might  be  due  to  an  underwater 
current,  to  surface  waves  moving  with  the  wind,  or  to  random  scatterer  motion. 

Spreading  loss,  also,  must  be  carefully  considered  in  the  modification  of  this 
reverberation  model.  In  SCATTER,  the  inverse  square  law  predicts  the  two-v/ay 
spherical-spreading  loss  as  being  proportional  to  R~*,  where  R  is  the  distance 
to  the  scatterer.  This  spreading  loss  is  constant  for  all  grid  cells.  When  a  non¬ 
isovelocity  environment  is  present,  some  grid  cells  may  encounter  a  focusing  of 
sound  and  others  a  diverging  of  sound.  An  alternate  method  for  determining  the 
spreading  loss  at  each  cell  has  been  used. 

A  series  of  raytracing  subroutines  has  been  written  to  perform  the  various 
tasks  necessary  to  construct  a  new  grid  over  a  constant  time  surfau:e.  These 
subroutines  also  compute  the  information  needed  to  determine  the  backscattering 
cross-section,  Doppler  shift,  frequency  spreading,  spreading  loss  and  absorption 
loss  for  each  cell. 

Whereas  REVMOD  does  not  allow  for  surface  or  bottom  reflections,  the 
revised  model  does.  Each  ray  is  traced  for  a  predetermined  length  of  time.  The 
rays  leading  to  a  particular  grid  cell  may  undergo  any  number  of  surface  and 
bottom  reflections.  These  reflections  are  assumed  to  be  specular  in  nature.  The 
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numbers  of  surface  and  bottom  reflections  are  determined  for  each  grid  ceil. 
Angularly  independent  loss  factors  for  the  surface  and  bottom  have  been  made 
input  variables  to  the  new  program  so  that  each  cell  may  be  weighted  by  a 
reflection  loss  parameter  not  previously  used  in  REVMOD. 

The  input  section  of  the  REVMOD  series,  INPUT,  has  undergone  minor 
changes  to  permit  sound  velocity  profile  and  other  pertinent  raytracing  input 
data.  Once  the  new  constant  time  grid  has  been  established  and  the  scattering 
function  determined  in  SCATTER,  the  remainder  of  the  REVMOD  series  of 
programs  operates  as  always,  except  that  the  modified  scattering  function  is 
incorporated. 


Chapter  2 

UNDERWATER  ACOUSTICS 
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2.1.  Reverberation 

2.1.1.  Introduction 

An  underwater  acoustical  environment  is  very  much  different  than  an  ideal 
free-space  acoustical  environment.  When  sound  is  transmitted  into  the  ocean 
from  an  underwater  source,  the  energy  encounters  various  inhomogeneities  in 
the  medium.  In  the  volume  of  the  ocean,  these  appear  as  marine  life,  such  as 
zooplankton,  passively  floating  animal  life,  and  fish.  Near  the  surface,  the  air- 
water  interface,  there  are  tiny  air  bubbles  and  a  rough  boundary.  The  bottom 
of  the  sea  is  not  smooth  either,  but  rather  is  made  up  of  peaks  and  valleys. 
These  discontinuities  cause  portions  of  the  transmitted  energy  to  be  scattered 
in  many  directions.  Some  of  this  scattered  energy  makes  its  way  back  to  the 
receiver,  which  is  at  the  source  location,  and  is  observed  as  reverberation.  The 
reverberation  level  at  any  given  time  after  the  so\md  has  been  transmitted 
is  dependent  upon  many  factors  in  addition  to  those  mentioned  previously. 
Included  are  the  transmission  loss  due  to  spreading  and  absorption,  reflection 
loss,  the  directivity  of  both  the  source  and  the  receiver,  and  the  abilities  of 
the  surface,  volume,  and  bottom  to  scatter  sound  back  to  the  receiver.  This 
scattering  of  sound  back  in  the  direction  from  which  the  sound  was  incident 
upon  a  scatterer  is  referred  to  as  backscattering. 

2.1.2.  Backscattering  Cross  Section 

A  measure  of  a  body’s  ability  to  scatter  sound  along  the  same  path  back  to 
the  receiver  at  the  source  location  is  called  the  backscattering  cross  section.  The 
backscattering  cross  section,  is  directly  related  to  the  cross-sectional  area 
(for  either  the  sxirface  or  bottom),  or  the  volume  of  the  scatterer.  A,  as  well  as 
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the  backscattering  function,  C6«  [2]- 

^6*  =  (2"1) 

The  backscattering  function  relates  the  backscattered  intensity  to  the  intensity 
incident  upon  a  scatterer.  The  scatterers  of  concern  are  either  surface,  volume,  or 
bottom  cells.  In  an  isovelocity  environment,  the  areas  and  volumes  of  these  cells 
are  easily  computed  geometrically.  The  zunount  of  symmetry  involved  leads  to  a 
small  number  of  calculations  required  to  define  the  areas  or  volumes  of  all  grid 
cells.  When  raytracing  is  involved,  however,  the  calculation  of  the  back-scattering 
area  or  volume  becomes  much  more  complicated  geometrically.  Although,  in  any 
one  horizontal  row,  the  size  of  the  grid  cells  is  constant,  the  grid  cells  may  vary 
in  size  within  a  vertical  column  for  non-isovelocity  conditions.  Thus,  many  more 
calculations  are  involved. 

For  the  determination  of  the  backscattering  volume  of  ocean  for  a  cell  at  a 
particular  distance,  a  pair  of  rays,  O  —  A  and  O  ~  B  in  Figure  2.1,  defining  the 
top  and  bottom  of  the  cell,  are  traced  for  a  given  amount  of  time,  t.  An  identical 
pair,  O  —  C  and  O  —  D,  are  traced  for  a  time  corresponding  to  the  back  of  the 
scattering  volume,  t  +  At,  where  At  is  to  be  determined  with  knowledge  of  the 
pulse  length.  The  endpoints  of  these  four  rays  define  an  area  in  a  vertical  plane, 
A  —  C  —  D  —  B.  The  scattering  volume  is  found  by  multiplying  this  area  by 
the  product  of  an  azimuthal  angular  increment,  A0,  and  a  radius,  r,  equal  to 
ct  cos  <^,  to  the  cell,  where  0  is  the  elevation  amgle.  See  Figure  2.1  for  an  example 
of  a  backscattering  volume  cell. 

A  special  case  arises  when  a  volume  cell  is  at  the  surface  or  bottom  of  the 
ocean.  In  the  general  case,  shown  in  Figure  2.2,  part  of  the  cell  has  just  undergone 
a  reflection  from  a  boundary  and  the  rest  of  the  cell  is  just  approaching  that 
boundary.  For  this  case,  the  cell  is  split  in  two  because  the  portion  of  the  cell 
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Figure  2.2  Volume  element  split  at  a  boundary. 
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after  the  reflection  would  have  undergone  a  boundary  forward  scattering  loss 
not  realized  by  the  first  portion  of  the  cell.  This  splitting  of  the  cell  is  done 
by  locating  the  centerline  of  the  ceil  and  dividing  it  where  it  impinges  upon 
the  boundary.  Each  portion  of  the  volume  cell  is  processed  separately  for  its 
contribution  to  the  reverberation,  amd  weighted  according  to  its  fraction  of  the 
total  cell  volume. 

Often,  the  predominant  reverberant  returns  come  from  the  ocean’s  surface 
and  bottom  as  opposed  to  the  volume.  Therefore,  it  is  very  important  to  be 
able  to  determine  the  areas  of  the  stirface  and  bottom  that  are  ensonified  at  a 
particular  time.  At  any  one  time,  there  may  be  many  different  patches  of  surface 
and  bottom  being  ensonified  due  to  different  raypaths  with  different  numbers 
of  reflections.  This  is  illustrated  in  Figure  2.3  by  two  constant  time  contours, 
one  at  the  front  and  or  e  at  the  back  of  the  scattering  volume.  The  contours 
are  defined  over  a  vertical  span  of  A^.  The  heavy  lines  at  the  surface  and  the 
bottom  represent  regions  of  boundary  ensonification.  One  area  of  the  surface  has 
been  ensonified  after  the  occurrence  of  a  bottom  reflection. 

2.1.3.  Scatterer  Depth 

The  backscattering  cross  section,  as  in  Equation  2.1,  cannot  be  determined 
unless  the  area  (in  the  case  of  either  the  surface  or  the  bottom)  or  the  volume 
of  the  scatterer  is  known.  For  a  volume  element  of  the  ocean,  the  A  in  Equation 
2.1  actually  represents  a  volume,  and  has  units  of  dBIm^.  Otherwise,  for 
ocean  surface  or  bottom  elements,  the  A  represents  an  area.  Whenever  a  pulse 
of  sound  b  transmitted,  the  ensonified  portion  of  the  medium  appears  as  a  shell 
with  a  thickness,  cr,  where  c  is  the  speed  of  sound  and  r  is  the  duration  of  the 
pulse.  In  the  case  of  an  isovelocity  medium,  the  ensonified  region  assumes  the 
geometry  of  a  spherical  shell.  The  thickness  corresponds  to  the  ph}rsical  length  of 
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Constant  time  contours  at  the  front  and  back  of  the 
scattering  volume  over  the  transmitted  span  of  angles, 
A^. 
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the  pulse  in  the  medium.  The  depth  used  for  the  thickness  of  a  scattering  element 
in  REVMOD,  however,  is  not  the  same  as  the  actual  ensonified  thickness  as  the 
pulse  travels  outward.  It  turns  out  to  be  one  half  of  this  thickness  as  discussed 
below. 

The  reverberation  observed  at  the  receiver  at  a  particular  time,  t,  after  the 
trailing  edge  of  a  pulse  has  been  transmitted,  is  the  net  result  of  all  energy 
reaching  the  receiver  at  that  time.  The  trailing  edge  of  a  pulse  is  transmitted 
one  pulse  length,  r,  later  than  the  leading  edge  of  the  pulse.  In  order  for  the 
backscattered  energy  from  the  leading  edge  of  the  pulse  to  arrive  at  the  receiver 
at  the  same  time  as  the  backscattered  energy  from  the  trailing  edge  of  the  pulse, 
the  leading  edge  must  travel  for  a  longer  period  of  time.  The  time  at  which 
the  trailing  edge  of  the  pulse  is  being  backscattered  is  which  is  half  of  the 
round  trip  time.  Since,  at  the  instant  the  trailing  edge  is  being  transmitted,  the 
leading  edge  has  already  been  traveling  for  one  pulse  length,  thb  leading  edge 
has  an  additional  “r”  added  to  the  total  travel  time,  t  as  is  shown  in  Figure 
2.4.  Therefore,  the  point  at  which  the  leading  edge,  line  B-B'  in  the  figure,  is 
backscattered,  such  that  it  reaches  the  receiver  at  the  saune  time  as  the  trailing 
edge,  line  A- A',  is  The  difference  in  these  travel  times  is 

f  +  r  t  _  ^ 

~~2  2  ~  2' 

'  Any  energy  received  as  reverberation  at  time  t  may  have  been  backscattered  at 

\  any  time  between  |  and  These  limits  define  the  scattering  depth,  or  shell 

thickness,  in  time  as  j. 

In  an  isovelocity  environment,  the  actual  depth  of  a  scattering  element,  in 
units  of  distance,  is  where  e  is  the  constant  sound  speed  in  the  medium.  When 
the  speed  of  sound  is  not  a  constant,  the  scatterer  depth  varies  throughout  the 
medium.  The  length  of  the  pulse,  in  time,  however,  remains  constant.  Thus,  it 
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Figure  2.4  Illustration  of  the  scatterer  depth  for  a  receiver  obser¬ 
vation  time,  t. 
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is  necessary,  when  non-isovelocity  raytracing  is  involved,  for  the  scattering  depth 
to  be  visualized  as  a  time  difference,  rather  than  a  physical  distance. 

2.1.4.  Mathematical  Description  of  Reverberation 

Underwater  reverberation  may  be  assumed,  as  a  first  approximation,  to  be 
the  result  of  acoustic  backscatter  from  various  discontinuities  in  the  oceanic 
environment.  Hodgkiss  [l]  has  described  the  ocean  as  a  linear,  r2mdom,  time- 
varying  filter,  through  which  a  transmitted  waveform,  /(t),  must  pass  before  it 
is  received  at  the  source  location  as  reverberation.  Random  variables  are  used 
because  there  is  seldom  enough  detailed  information  on  the  acoustic  scattering 
environment  to  use  a  deterministic  description.  In  this  section,  the  concept  of 
a  scattering  function  is  developed.  The  scattering  function  describes  how  the 
energy  from  a  transmitted  pulse  is  statbtically  redistributed  over  frequency  and 
time  [3]. 

Van  Trees  [4]  has  shown  how  a  finite  length  scatterer  may  be  modeled  as  a 

multi-element  reflector.  The  scatterer,  of  length  L  in  time,  is  broken  down  into 

increments  of  length  AA  as  is  shown  in  Figure  2.5.  The  total  return  from  this 

scatterer  may  be  expressed  as  a  sununation  of  returns  from  all  increments. 

N 

i(t) = (2-3) 

»=0 

where  i(t)  b  the  backscattered  waveform,  Et  b  the  source  energy  level,  f[t)  is  the 
transmitted  waveform,  and  b{Xi)  b  a  complex  Gaussian  weighting  variable.  As 
the  size  of  the  increment  tends  toward  zero,  the  summation  becomes  an  integral 


w. 
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=  fit-\)bR{x)dx 

Lo 


(2.4) 


In  Equation  2.4,  ^^(A}  b  a  sample  function  from  a  zero-mean  complex  Gaussian 
process. 
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Figure  2.5  Geometry  of  a  scattering  element. 
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Hodgkiss  has  used  this  theory  to  model  reverberation  as  a  result  of  the 
acoustic  backscatter,  i(t),  from  the  ocean.  He  has  implemented  a  random,  time- 
varying,  complex,  impulse  response  of  a  filter,  6(t  -  j,A),  and  has  obtained  the 
following  [l]: 

4-00 

5(0  =  \/^  f  j,X)dX.  (2.5) 

— oo 

Here,  A  is  a  time  delay  variable  which  corresponds  to  distance.  The  time  variation 
of  the  impulse  response  function  has  a  covariance  function  of  the  form 

^j(t-u,A).  (2.6) 

Since  b(t  —  ^,A)  has  been  assumed  to  be  a  zero-mean  process,  the  covariance 
function  is  equal  to  the  correlation  function  [5].  A  spectral  density  function  is 
obtained  through  the  Fourier  transform  of  the  correlation  function  [5].  Hodgkiss 
has  described  the  scattering  function  as  the  Fourier  transform  of  the  covariance 

function  associated  with  the  impulse  response  of  the  filter  [l], 

+00 

Si(/,X)=  I  (2,7) 

— OO 

The  scattering  function  describes  not  only  the  attenuation  due  to  propagation 
losses  and  source  and  receiver  directivities,  but  also  the  redistribution  in 
frequency  due  to  scatterer  motion. 

2.1.5.  Computer  Modeling  of  Reverberation 

In  section  2.1.4,  a  mathematical  description  of  reverberation  was  discussed 
in  which  the  ocean  was  approximated  as  a  complex,  time-varying  filter  through 
which  a  transmitted  signal  must  pass  in  order  to  be  received  as  reverberation. 
The  time-varying  impulse  response  of  the  filter,  b{t  —  ^,A),  would  be  difficult, 
if  not  impossible,  to  predict  in  the  computer  modeling  of  a  realistic  oceanic 
environment.  However,  the  concept  of  the  scattering  function,  5j(/,  A),  may 
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quite  easily  be  used  to  develop  a  practical  model  for  reverberation.  Hodgkiss 
has  shown  how  the  fundamental  laws  of  acoustics,  zis  well  as  some  geometric 
considerations,  caji  be  used  to  construct,  piece  by  piece,  a  scattering  function  of 
the  ocean. 

At  any  distance,  the  ocean  may  be  divided  into  a  number  of  sections,  or 
cells,  each  contributing  to  the  total  scattering  function  for  that  distance.  It  is 
convenient  to  divide  an  isovelocity  environment  into  rows  and  columns  falling 
on  the  surface  of  a  sphere  with  a  raditis  equal  to  the  scatterer  distance  from  the 
source.  Therefore,  a  particular  cell  may  be  referred  to  as  being  located  in  the 
row  and  the  colunm  for  each  distance.  An  incremental  scattering  function, 
for  each  cell,  in  which  both  Doppler  spreading  and  frequency  shifting  axe  taken 
into  account,  is  produced.  The  incremental  scattering  functions  each  contribute 
to  the  total  scattering  function. 

Doppler  shift  for  a  source  moving  with  some  velocity,  v,  at  an  angle,  7,  with 
respect  to  the  scatterer,  is  calculated  through  the  use  of  the  Equation  [6]: 

r  j-  /,  2i;cos7. 

fo  =  /r(i  H - -)  (2.8) 

c 

where  //j  is  the  Doppler  shifted  frequency,  /y  is  the  transmitted  frequency,  and  c 
is  the  speed  of  sound,  which  has  thus  far  been  assumed  to  be  constant  throughout 
the  ocean.  See  Appendix  A.l  for  a  derivation.  An  additional  Doppler  shift  arises 
due  to  a  cell  moving  with  a  mean  velocity,  u,  at  an  angle,  7,  with  respect  to  the 
source.  The  equation  governing  this  case  turns  out  to  be  identical  to  Equation 
(2.8),  except  that  v  in  Equation  (2.8)  is  replaced  by  u.  See  Appendix  A.2  for  a 
derivation.  The  net  Doppler-shifted  component  is  the  summation  of  the  Doppler 
due  to  a  moving  source  and  the  Doppler  due  to  a  moving  scatterer. 
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The  Doppler  shifting  discussed  to  this  point  pertains  only  to  the  monostatic 
case,  in  which  the  sound  is  received  at  the  same  location  from  which  it  originated. 
The  question  hzis  arisen  as  to  whether  a  specular  forward  reflection  from  a  moving 
scatterer  would  contribute  a  component  of  Doppler  shift.  Pierce  [7]  has  shown 
how  sound  scattered  in  any  direction,  not  necessarily  backward,  is  influenced  by 
the  Doppler  effect.  A  situation  in  which  the  sound  is  not  backscattered,  but  is 
scattered  in  some  other  direction,  may  be  referred  to  as  a  bistatic  case  because, 
in  general,  an  observation  point  from  which  the  Doppler  shift  is  realized,  is  not 
located  at  the  source.  The  geometry  for  a  case  of  bistatic  scatter  is  shown  in 
Figure  2.6.  Va  is  the  velocity  vector  of  the  scatterer,  and  n,-  and  n,  are  unit 
vectors  in  the  directions  of  the  incident  and  scattered  waves  respectively.  The 
vector,  ej,-,  is  the  unit  vector  in  the  direction  n,  —  n,-.  Pierce’s  equations  may  be 
rearranged  to  show  that 

/d  = /r(l- —•«6»2siniA'r).  (2.9) 

C  i 

In  Equation  2.9,  c  is  the  speed  of  sound  at  the  scatterer  and  A7  is  the  deflection 
angle  between  n,-  and  n^.  If  the  bistatic  reflection  were  to  be  in  the  form  of  a 
specular  reflection,  and  the  velocity  vector,  Vg,  were  assumed  to  lie  in  the  x-y 
plane,  then  the  vectors,  Vg  and  would  be  orthogonal  and  the  dot  product  of 
the  two  would  equal  zero.  Figure  2.7  illustrates  this  case.  Therefore,  in  the  case 
of  a  forward  speculeir  reflection  from  a  scatterer  moving  horizontally,  such  as  the 
surface  of  the  ocean,  no  frequency  shifting  would  result. 

Frequency  spreading  for  a  scattering  element  of  the  ocean  is  a  result  of 
variations  in  the  scatterer  speed  from  the  mean  value.  The  scatterer  speed 
distribution  has  been  assumed  to  be  a  zero-mean  Gaussian  probability  density 
function,  which  is  characterized  by  its  standard  deviation,  o.  The  Gaussiaui 
probability  density  function,  which  is  used  to  closely  approximate  many  random 
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phenomenA,  m^y  be  exprenc<l  m  iSl, 

p(r)  =  — (2.10) 
aV2» 

where  z,  the  KAtterer  speed,  U  the  raAdoro  variable,  a  is  the  si&ndard 
deviation,  and  p  ia  the  mean  scatterer  speed  which  has  been  assumed  equal 
to  zero  for  spectral  spreading  considerations.  The  random  spectral  spreading  is 
directly  related  to  the  random  scatterer  motion  through  the  Doppler  equation 
shown  previously.  Both  spectral  spreading  and  Doppler  shifting  comprise  the 
incremental  scattering  function  for  each  cell. 

The  previously  mentioned  redistribution  of  frequency  between  the  lime 
of  transmitting  and  receiving,  for  each  ocean  element,  is  referred  to  as  the 

S 

incremental  scattering  function.  (ij.  Incremental  scattering  functions 

must  be  calculated  for  each  cell  and  weighted  by  an  attenuation  factor, 
before  being  summed  to  obtain  the  composite  scattering  function.  The 
attenuation  for  each  scattering  element  is  comprised  of  absorption  loss,  spreading 
loss,  and  the  directivities  of  the  source  and  receiver.  The  backscattering  strength 
has  been  assumed  to  be  independent  of  the  angle  of  incidence.  Hodgkiss  has  not 
considered  either  siirface  or  bottom  forward  reflections  in  his  model,  thus,  no 
energy  is  assumed  to  be  lost  due  to  interactions  with  these  boundaries.  The 
backscattering  strength  of  each  cell  is  also  taken  into  account  in  the  attenuation 
factor  through  the  backscattering  coefficient.  The  total  scattering  function  for 
the  ocean,  then,  may  be  expressed  as  a  summation  of  all  incremental  scattering 
functions,  each  weighted  by  a  corresponding  attenuation  factor  [l]. 

«,(/.*)  A)  (2-11) 

The  total  backscattered  power  levels  may  thus  be  expressed  as  a  summation  of 
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the  attenuation  factors  [ij. 

^g(0,A)  =  ^A,-y  (2.12) 

Equation  2.12  is  used  to  determine  the  total  power  gain  verstis  time  for 
a  particular  scenario.  The  scattering  level  for  each  cell,  BSL^  j,  is  computed 
by  first  finding  the  volume  of  the  cell,  or  the  area  of  the  cell  if  it  falls  on  the 
surface  or  the  bottom,  and  then  multiplying  by  an  appropriate  backscattering 
coefficient,  in  units  of  either  dBfm^  or  dB/m},  depending  on  the  type  of  cell 
being  processed.  Absorption  loss  and  spreading  loss  factors  are  included  in  the 
attenuation  factor  equation  as  are  weighting  factors  corresponding  to  the  source 
and  receiver  directivities  [l].  The  attenuation  factors  may  be  expressed  as  follows. 

10  log  A,- y  =  BSLi^j  {dB)  -  ALij  {dB)  -  SL^^j  {dB)  -  BAij  {dB),  (2.13) 

where  BSLls  the  backscattering  level,  AL  and  SL  are  absorption  and  spreading 
loss,  respectively,  and  BA  is  the  combined  transmit  and  receive  beam  pattern 
weighting.  Equation  2.13  makes  up  part  of  the  SONAR  (SOund  NAvigation 
Ranging)  equation  [l].  The  complete  SONAR  equation  includes  additional  terms 
related  to  the  transmitter  and  receiver,  such  as  sensitivities,  transmitted  power 
level  and  the  ambient  noise  in  the  ocean.  These  factors  are  out  of  the  scope  of 
this  effort  and  are  irrelevant  to  the  scattering  function  determination. 

2.2.  Sound  Propagation:  Raytracing 

When  acoustic  waves  travel  from  one  fluid  to  another  with  a  different  speed 
of  sound,  c,  at  some  angle  other  thsm  normal,  the  waves  are  bent,  or  refracted. 
The  sound  propagates  through  the  second  medium  at  a  different  angle  than  it 
did  through  the  first.  If  a  number  of  fluids  with  varying  speeds  of  sound  are 
layered,  a  wave  will  refract  through  each  layer.  Tracing  a  wave  through  such  a 
varying  medium  can  be  made  much  easier  by  the  use  of  ray  theory  as  opposed  to 
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wave  theory.  Sound  rays  are  defined  as  lines  which  are  everywhere  perpendicular 
to  acoustic  wavefronts,  i.e.,  areas  of  constant  phase  [6].  Ray  theory  is  not  an 
exact  replacement  for  wave  theory,  but  is  a  good  approximation  under  certain 
conditions. 

2.2.1.  The  Eikonal  Equation 

Kinsler,  Frey,  Coppens,  and  Sanders  [6]  provide  a  derivation  of  the  Eikonal 
equation  by  assuming  a  solution  of  the  following  form  to  the  wave  equation: 

p{x,y,z,t)  =  (2.14) 

where  P  is  pressure,  F  is  distance,  and  Co  is  a  constant  value  of  phase  speed. 
The  speed  of  sound,  c,  is  assumed  to  be  a  function  of  x,  y,  and  z.  Areas  over 
which  r  =  constant  are  surfaces  of  constant  phase.  VF  is  always  perpendiculau^ 
to  these  constant  phase  surfaces  [6].  The  operator,  V,  is  defined  as 


„  d^.  d  *.  d  t 

+  al*  • 


Likewise,  the  operator,  V^,  is  defined  as 


5* 

5x*  5y*  dz^ 

When  the  trial  solution  is  substituted  into  the  wave  equation,  the  result  is 

^  ■  vr  +  v>r)  =  0.  (2.15) 


All  terms,  except  for  the  second  and  third,  can  be  considered  very  small  if  both 
P  and  vr  are  slowly  varying.  Neglecting  them  yields  the  Eikonal  equation: 


vr  ■  vr  i 


(2.16) 


where 


n(x,y,s) 


Co 

c(x,y,z) 


(2.17) 
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is  the  refractive  index  and  Co  is  an  arbitrary  reference  phase  speed.  There  are 
two  conditions  necessary  for  Equation  2.15  to  reduce  to  Equation  2.16.  The  first 
is  that  the  amplitude  of  the  wave  must  not  change  appreciably  over  distances 
comparable  to  a  wavelength.  The  second  is  that  the  speed  of  sound  must  not 
change  appreciably  over  distances  comparable  to  a  wavelength.  The  solution  of 
the  Eikonal  equation  for  F  yields  the  path  that  the  acoustic  energy  follows. 

2.2.2.  Snell’s  Law 

As  a  ray  passes  from  one  medium  into  another,  the  angle  of  refraction  is 
determined  by  both  the  angle  of  incidence  and  the  speeds  of  sound  of  the  two 
layers.  This  relation  is  given  by  Snell’s  Law  (2): 

sin^i  _  sin  ^2 

Cl  C2 

In  Equation  2.18,  ci  is  the  speed  of  sound  in  the  first  medium,  and  C2  is  the 
speed  of  sound  in  the  second  medium.  Likewise,  and  ^2  the  angles  of  the 
rays  in  the  first  and  second  media,  respectively.  See  Figure  2.8  for  an  illustration 
of  a  refracting  ray  and  the  corresponding  wavefronts.  For  the  case  shown,  ci  > 
C2.  Snell’s  Law  also  applies  in  a  multi-layered  medium.  As  the  layer  thickness 
becomes  infinitesimal,  the  speed  of  sound  becomes  a  continuous  function  of  depth. 
The  angle,  <f>,  at  a  depth,  a,  of  a  ray  beginning  at  an  angle,  4>i,  and  depth,  z,, 
cam  be  determined  by 


(2.18) 


sin<^  _  sin<^» 
'(»)  ~  '(*i) 


(2.19) 


where  a  is  the  Snell’s  Law  constant  of  the  ray.  The  angles  are  always  measured 
between  the  ray  and  the  axis  normal  to  the  layers. 

Clay  and  Medwin  [2]  show  that  if  the  medium  is  divided  into  horizontal 
stratifications  of  constant  velocity  gradient,  dcfdz,  then  the  following  expressions 
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Figtire  2.8  Refracting  ray  and  wavefronts.  (ci  >  cj). 
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for  time  and  distance  may  be  obtained: 


1.  _ 


(2.20) 


Where 


rf-ri  =  -j(cos*i  -cos*^) 


(2.21) 


ti  =  initial  time  of  a  ray  in  a  layer, 

tf  =  final  time  of  the  ray  in  the  sarnie  layer, 

Ti  =  initial  horizontal  distance  of  a  ray  in  a  layer, 
rf  =  final  horizontal  distance  of  the  ray  in  the  same  layer, 
b  =  sound  velocity  gradient  in  layer, 
a  =  Snell’s  Law  constant  for  the  ray, 

<f>i  =  initial  angle, 

<f>f  =  final  amgle, 

Wi  =  transformed  initial  depth  variable, 

Wf  =  transformed  final  depth  variable. 

The  transformed  depths  are  as  follows: 

Wi  =  Zi  -zi  +  (2.22) 

Wf  =  Zf  -  Zi  +  ^^^  (2.23) 

where  c{z\)  is  an  arbitrary  reference  speed  of  sound  at  depth  zi.  These  equations 
are  valid  for  travel  times  and  distances  through  one  particular  layer.  A  ray  must, 
therefore,  be  traced  through  one  layer  at  a  time.  The  rays  which  travel  through 
constant  velocity  gradient  layers  follow  the  arcs  of  circles  [8].  The  radius  of  the 
circle  is  equal  to  {ab)~^  and  is  constant  for  a  particular  ray  in  a  particular  layer. 
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2.2.3.  Sound  Velocity  Profile 

Due  to  the  tendency  of  less  dense  water  to  rise  above  more  dense  water,  the 
ocean  may  be  modeled  as  a  horizontally  stratified  medium.  Therefore,  the  speed 
of  sound  may  be  considered  to  vary  only  with  depth.  This  might  not  be  true  in 
some  situations,  especially  over  long  distances.  However,  this  is  the  assumption 
used  in  the  development  of  the  raytracing  algorithms  for  the  modification  of  the 
isovelocity  reverberation  model.  A  plot  of  the  speed  of  sound  versus  depth  is 
referred  to  as  a  sound  velocity  profile. 

The  speed  of  sound  is  usually  measured  by  the  use  of  either  a  velocimeter 
or  a  bathythermograph.  A  velocimeter  is  an  acoustical  device  for  measuring  the 
actual  sound  speed.  Urick  [9]  gives  a  brief  description  of  Tschiegg  and  Hays’ 
[10]  original  velocimeter.  The  bathythermograph  measures  the  temperature  as  a 
function  of  depth  [9].  With  some  knowledge  of  or  an  assumption  of  the  sadinity, 
and  a  calculation  of  the  pressure,  the  temperature  may  be  used  to  determine  the 
speed  of  sound  through  Wilson’s  equation  [11]. 

V  =  1449.22  +  AVt-  +  AVp  +  AV^  +  AV^j-p  (2.24) 

where 

AVt  =  4.572ir  -  4.4532  X  lO'^r*  -  2.604  x  10~‘‘r®  +  7.9851  x  10“®r'‘, 

AVp  =  1.60272  X  10-^P  +  1.0268  x  lO'^P*  +  3.5216  x  lO'^P®, 

-3.3603  X 

AVs  =  1.39799(5  -  35)  +  1.69202  x  10-3(5  -  35)3, 

and 

AVstp  =  {S-  35)(-1.1244  x  lO'^T  +  7.7711  x  lO-^T*  +  7.7016  x  IQ-^P 
-1.2943  X  10-3p2  +  3.1580  x  lO'^PT  +  1.5790  x  IQ-^PT^) 
+P(-1.8607  X  10-‘‘r  +  7.4812  X  IQ-^T^  +  4.5283  x  10-37®) 
+p3(-2.5294  X  10-’'r+ 1.8563  x  lO'^T®)  +  p3(- 1.9646  X  10- ^°r). 
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A  similar  formulation  was  made  earlier  by  Wilson  [12]  for  a  limited  range  of 
salinity  (33ppt  <  S  <  37ppt).  Here,  ppt  stands  for  parts  per  thousand.  Wilson 
had  later  revised  that  formulation  to  obtain  Equation  2.24,  which  is  valid  in  the 
range  of  0  <  5  <  37ppt.  Equation  2.24  is  also  valid  for  the  temperature  range  of 
— 4°C  <T  <  30°C  and  the  pressure  range  of  1kg  f  cm?  <  P  <  lOOOkgfcm^  [ll]. 
Actually,  the  units  given  for  pressure  are  incorrect.  A  factor  for  the  acceleration 
due  to  gravity,  g,  is  missing.  If  g  is  assumed  to  be  equal  to  9.81m/s^,  then 
the  range  of  pressures  for  which  Equation  2.24  is  valid  is  9.81  x  lO^Pa  <  P  < 
9.81  X  lO^Pa.  The  temperature,  pressure,  and  salinity  ranges  of  validity  may  be 
expanded  slightly,  introducing  only  a  small  amount  of  error  [13|. 

It  can  be  seen,  then,  that  the  speed  of  sound  increases  with  an  increase  in 
all  three  parameters  -  temperature,  pressure,  and  salinity.  The  sound  velocity 
profile,  th\is,  changes  with  temperature,  pressure,  and  salinity.  At  any  given 
location,  the  sound  velocity  profile  varies  with  time  both  throughout  the  day 
and  with  the  seasons  [9]. 

2.2.4.  Ray  Diagrams 

Raytracing  plots,  or  ray  diagrams,  may  be  used  as  helpful  visual  aids 
in  understanding  the  behavior  of  acoustic  fields.  A  ray  diagram  is  a  set  of 
rays,  usually  angularly  equispaced  at  the  source,  which  dse  plotted  with  depth 
as  a  function  of  the  horizontally  traveled  distance.  Differences  in  intensity 
levels  at  different  points  in  the  medium  may  be  observed  by  comparing  the 
spacing  of  adjacent  rays  at  some  distance  away  from  the  source.  Shadow  zones, 
regions  where  there  is  very  little  acoustic  intensity,  and  caustics,  regions  of  very 
concentrated  intensity,  may  easily  be  located  in  raytracing  plots.  Figure  ^.9 
shows  a  ray  diagram  which  is  characteristic  of  an  isovelocity  environment.  The 
rays  have  been  plotted  at  1.0®  increments  for  the  range  of  —5.0®  <  4>  <  +5.0®, 
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where  4>  is  the  vertical  angle,  in  degrees,  of  the  rays  at  the  source.  The  angle,  4>, 
is  measured  from  the  horizontal  plane  in  which  the  source  lies.  Positive  angles 
are  assumed  in  the  upward  direction  from  this  plane  and  negative  angles  are 
assumed  in  the  downward  direction.  All  rays  are  straight  lines,  which  should 
be  expected,  as  no  refraction  occurs  in  an  isovelocity  environment.  Because  the 
ray  spacing  is  constant  at  any  particulau'  distance,  the  intensity  is  uniform  at  all 
points  in  the  medium  at  any  fixed  distance  from  the  source. 

In  Figure  2.10,  a  non-isovelocity  sound  velocity  profile  was  assumed.  This 
profile  is  shown  in  Figure  2.11.  The  rays  that  have  been  traced  are  the  same 
as  those  that  appear  in  Figure  2.9.  It  is  evident  in  this  diagram  that,  in  any 
particular  region,  the  spacing  between  the  rays  is  not  constant.  This  is  an 
indication  that  the  intensity  is  not  constant  for  ail  points  at  the  same  travel 
time  from  the  source.  In  region  A,  the  rays  are  closely  spaced,  indicating  a 
relatively  high  intensity,  whereas,  in  region  B,  the  rays  are  separated,  indicating 
a  lower  intensity.  It  can  also  be  seen  in  Figure  2.10  that  some  of  the  sound  rays 
beginning  upward  from  the  source,  refract  away  from  the  surface,  as  opposed  to 
the  ones  that  refiect  from  the  surface  in  Figure  2.9. 

An  alternate  way  of  examining  the  behavior  of  sound  through  raytracing  is 
to  construct  the  actual  wavefronts  at  different  points  in  time  after  the  sound 
has  been  transmitted.  This  can  be  done  by  tracing  a  series  of  rays,  each  for 
the  same  amount  of  time,  and  retaining  only  the  final  point  of  each  ray.  These 
points  represent  one  instant  in  time,  thus,  they  may  be  connected  to  form  a 
wavefront.  A  series  of  these  wavefronts  may  be  plotted  to  illustrate  how  the 
acoustic  waves  travel  over  time.  Shadow  zones  and  the  general  refraction  of 
sound  may  be  observed  rather  easily  with  this  type  of  plot,  although,  there  is  no 
major  advantage  over  the  ray  diagrams.  Figure  2.12  shows  a  series  of  wavefronts 
superimposed  on  the  associated  ray  diagram.  An  isovelocity  scenario  was  used 
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Figure  2.10 


Ray  diagram  illustrating  regions  of  varying  intensity 
for  a  non-isoveiocity  case. 
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Figure  2.12  Wavefronts  superimposed  on  a  ray  diagram  for  the 
isovelocity  case. 
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and  the  rays  were  traced  from  -5.0°  <  ^  <  -t-S.Oo  at  1.0°  increments.  It  can 
be  seen  that  the  wave  has  undergone  a  surface  reflection,  after  which  a  reflected 
wavefront  has  appeared.  Because  of  the  different  scaling  on  the  z  and  y  axes, 
the  wavefronts  do  not  appear  to  be  curved  as  they  are. 

In  Figure  2.13,  a  decreasing  speed  of  sound  with  depth  has  been  assumed, 
and  the  z  axis  has  been  contracted.  The  wavefronts  have  been  plotted  from 
-10.0°  <  <  +10.0°,  as  well  as  the  ray,  at  approximately  8.0°,  that  grazes 

the  surface.  The  formation  of  a  surface  reflection  and  a  shadow  zone  are  both 
illustrated  in  this  figure. 

2.2.5.  Ray  State  Numbers 

Due  to  Snell’s  law  and  the  assumption  of  a  horizontally  stratified  medium, 
any  ray  beginning  at  an  angle  above  the  horizontal  will  eventually  return  to 
the  source  depth  at  the  same  angle  below  the  horizontal.  The  ray  may  either 
reflect  specularly  from  the  surface  of  the  ocean  or  it  may  refract  downward, 
i.e.,  turn  before  coming  into  contact  with  the  surface.  The  same  principle  of 
symmetry  applies  for  rays  beginning  downward  from  the  source  which  either 
reflect  or  refract  upward,  back  to  the  source  depth.  An  ocean  bottom  with  a 
constant  depth  has  been  assumed.  See  Figure  2.14  for  examples  of  raytracing 
symmetry.  All  ray  segments  labeled  A,  have  identical  travel  times  and  distances. 
The  same  applies  to  the  travel  times  and  distances  for  the  segments  labeled,  B. 

As  an  aid  to  the  discussion  of  different  ray  paths  from  a  source  to  a  scatterer, 
the  concept  of  ray  state  numbers  will  be  introduced.  Let  the  ocean  be  divided  into 
three  layers,  the  first  being  the  region  between  the  bottom  and  the  source;  the 
second,  that  between  the  source  and  the  scatterer;  and  the  third,  that  between 
the  scatterer  and  the  ocean  surface.  It  may  be  assumed  that  the  source  is 
always  below  the  scatterer  for  purposes  of  labeling  individual  raypaths.  Acoustic 
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Figure  2.13  A  negative  velocity  gradient  produces  downwardly 
refracting  waves.  A  shadow  zone  is  formed  at  the 
stirface. 
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Figtxre  2.14  Raytracing  symmetry. 
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reciprocity  dictates  that  the  roles  of  the  source  and  the  receiver,  i.e.,  scatterer, 
may  be  interchanged. 

Roeckel  [I4j  has  developed  a  method  of  ordering  rays  according  to  the  number 
of  times  they  pass  through  the  region  between  the  source  and  the  scatterer.  This 
ntunber  of  passes  is  denoted  by  the  state  number,  KD.  Since  all  rays  begin  at 
the  source  and  end  at  the  scatterer,  KD  must  always  be  an  odd  number.  Each 
class  of  ray  solutions  is  represented  by  the  index,  K,  such  that 

KD  =  2K-  1.  (2.25) 

Two  other  state  numbers,  KL  and  KA,  have  also  been  defined.  The  number  of 
times  a  ray  passes  through  the  region  below  the  source  is  represented  by  KL. 
Likewise,  the  number  of  times  a  ray  passes  through  the  region  above  the  source 
is  represented  hy  K A.  KL  and  KA  are  each  related  to  K  in  one  of  two  ways, 
depending  on  the  actual  raypath. 

KL  =  2K  or  2(K  ~  1).  (2.26) 

KA  =  2K  or  2{K  -  1).  (2.27) 

See  Appendix  B  for  the  class  of  rays  for  iC  =  1.  For  any  index,  K,  a  particular 
ray  is  uniquely  defined  by  a  set  of  state  numbers,  KD,  KL,  and  KA. 

2.2.6.  Final  Position  and  Angle  of  Ray 

The  final  position  and  the  final  angle  of  a  ray  from  the  source  after  a  certain 
length  of  time  are  important  parameters  to  be  determined.  The  final  position 
of  a  ray  will  allow  for  the  calculation  of  the  grid  cell  thickness,  as  well  as  the 
ensonified  volume  or  area  of  any  grid  cell.  The  final  angle  of  the  ray  at  the 
scattering  element  will  determine  the  amount  of  Doppler  shift  contributed  by 
that  element. 
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The  procedure  used  to  calculate  the  final  position  of  a  ray,  given  a  total  travel 
time,  t,  consists  of  two  basic  steps.  First,  the  layer  in  which  the  ray  ends  after 
time,  t,  has  passed  is  found.  Then,  the  actual  position  of  the  ray  within  this 
layer  is  determined.  In  order  for  this  to  be  done,  an  equation  which  yields  a  final 
depth  given  an  initial  angle,  depth  and  time  must  be  developed.  This  can  be 
achieved  by  beginning  with  Clay  and  Medwin’s  [2]  alternate  expression  for  the 
time  traveled  in  a  layer: 


_  _1  ii;/(l-h(l 

^  +  (1  - 


(2.28) 


(2.29) 


Rearranging,  and  solving  for  Wf,  the  transformed  depth  variable,  the  following 
result  is  obtained: 

r  1 

Wf  =  2 - ^ — - 5 — -  . 

The  final  depth  of  the  ray  in  the  layer  is  obtained  by  transforming  Wf  back  to 
the  actual  depth,  Zf. 

Zf  =  Wf  +  zi  —  — (2.30) 


Since  these  equations  are  valid  for  only  one  layer  at  a  time,  it  is  necessary  to  find 
the  layer  in  which  the  solution  lies  before  calculating  the  exact  final  position  of  the 
ray.  Once  a  final  depth,  Zf,  is  obtained,  the  final  angle  and  horizontally  traveled 
distance  follow  through  the  use  of  Snell’s  law  and  Equation  2.21,  respectively. 

The  total  distance  traveled  along  the  ray  is  needed  for  the  determination  of 
the  absorption  loss.  Since  linear  velocity  gradient  segments  are  usckl  in  the  sound 
velocity  profile,  the  path  of  any  ray  through  any  layer  will  always  follow  the  arc 
of  a  circle,  as  was  mentioned  earlier.  The  arclength  is  simply  the  angle,  V’l  times 
the  radius,  (a6)~^.  See  Figure  2.15  for  the  geometry  used  in  determining  the 
arclength  of  a  ray  passing  through  any  one  layer.  From  Figure  2.15,  it  can  be 
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seen  that 


thus, 


IT 

02  =  --  <^2, 


tj}  =  01  —  02  =  <i>2  —  <i>i. 


and  finally, 

ARCLENGTHi  =  (2-31) 

ao^ 

where  the  subscript  i  denotes  layer  i  and  all  angles  are  in  radians.  As  a  ray  is 
traced  through  various  layers  along  its  path  from  source  to  scatterer,  the  total 
horizontal  distance  and  arclength  may  be  incrementally  obtained  through  the 
use  of  Equations  2.21  and  2.31  for  each  layer. 


2.2.7.  Critical  Angles 


A  critical  angle  is  defined  here  as  an  angle  at  the  source,  such  that  the 
corresponding  ray  reaches  a  turning  point  at  a  depth  at  which  the  gradient  of 
the  speed  of  sound  changes.  This  is  any  depth  dividing  two  constant  velocity 
gradient  bands  in  the  ocean.  A  turning  point  is  defined  as  a  depth  at  which 
a  ray  changes  vertical  direction  from  downward  to  upward  or  from  upward  to 
downward.  These  critical  angles  are  essential  for  classifying  different  types  of 
ray  paths.  All  rays  v/hich  fall  between  two  adjacent  critical  angles  will  behave 
similarly  and  will  travel  together  in  a  group. 

The  only  depth  at  which  a  ray  will,  or  can,  turn  is  one  at  which  the  speed  of 
sound  is  greater  than  or  equal  to  the  speed  of  sound  at  the  source.  This  can-  be 
demonstrated  by  the  use  of  Snell’s  law  and  the  fact  that  when  a  ray  turns,  its 
direction  is  horizontal,  or  ninety  degrees  with  a  reference  of  zero  degrees  being 
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straight  downward: 

sin<^c  _  sinOO  _  1 

Csouree  Clurn  C(urn 

Therefore, 


(2.32) 


(2.33) 


Thus,  a  critical  angle,  <i>c,  can  only  exist  when  the  speed  of  sound  at  the  turning 
point  is  greater  than  or  equal  to  that  at  the  source. 

Since  a  ray  beginning  at  some  angle  below  the  source  will  eventually  be 
directed  at  an  equal  angle  above  the  source  and  vice  versa,  every  critical  angle 
has  a  mirror-image  angle,  which  is  also  critical,  associated  with  it.  See  Figure 
2.16  for  examples  of  rays  at  critical  angles.  Note  that  all  points  on  the  sound 
velocity  profile  with  a  speed  of  sound  greater  than  the  speed  of  sound  at  the 
source  depth  will  yield  a  critical  angle  both  above  and  below  horizontal. 


2.2.8.  Losses 


Along  its  path  from  the  source  to  the  scatterer  and  back  to  the  receiver, 
the  transmitted  energy  encounters  three  types  of  attenuation,  or  transmission 
loss  {TL)\  spreading  loss,  absorption  loss,  and  reflection  loss.  Non-isovelocity 
propagation  can  have  a  great  influence  on  these  loss  mechanisms.  In  a  non- 
isovelocity  environment,  not  only  does  the  absorption  coefficient  vary  with  depth, 
but  the  spreading  loss,  to  the  scatterer  and  back,  is  no  longer  R~*.  Reflection 
losses  due  to  interaction  with  the  surface  and  the  bottom  aure  also  affected  by 
ray  bending.  Under  certain  conditions,  the  sound  velocity  profile  may  cause  the 
acoustic  energy  to  be  refracted  away  from  the  surface  or  bottom  in  which  case 
no  reflection  loss  would  occur. 

The  acoustic  intensity  at  a  field  point  in  a  non-isovelocity  environment  may 
be  approximated  through  the  use  of  ray  theory.  Depending  on  the  sound  velocity 
profile,  the  source  location  and  the  field  point,  the  rays  may  either  converge  or 
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diverge.  Clay  and  Medwin  [2j  show  that  the  separation  distance  of  closely  spaced 
rays  in  the  vicinity  of  an  observation  point  caoi  be  used  to  estimate  the  pressure 
at  that  point  from  an  omnidirectional  source.  This  pressure  may  be  expressed  in 


the  form: 


where, 


rL  PoCo 


p2  _  PqCqTI 

^  ~  47r/22  • 


(2.34) 


See  Figxire  2.17  for  an  illustration  of  the  geometry.  In  Equation  2.34,  ^  =  0 
corresponds  to  an  angle  straight  downward,  and  A(^  must  be  in  radians.  The 
transmission  loss  due  to  spreading  is  defined  [2j  as 


(2.35) 


In  terms  of  acoustic  intensity,  this  may  be  represented  as 


(2.36) 


The  expression  for  the  squared  pressure  at  a  field  point  (Equation  2.34) 
includes  a  term  which  is  the  ratio  of  the  density  at  the  field  point  to  that  at  the 
source.  Since  the  density  changes  just  a  fraction  of  a  percent  with  significant 
changes  in  temperature,  pressure,  and  salinity  [15],  the  minute  change  in  density 
can  be  neglected  in  the  pressure  calculations.  Clay  and  Medwin  [2]  point  out 
that  L  in  Figure  2.17  may  be  estimated  from  the  angle,  at  the  receiver  and 
numerical  values  of  h  if  the  rays  are  not  converging  or  diverging  rapidly.  This 


relation  is  given  as 


L  =  jhsin^j. 


(2.37) 


If  two  closely  spaced  rays  are  traced  for  the  same  amount  of  time,  the  straight-line 
separation  of  the  rays  may  be  solved  exactly.  Since  the  two  rays  are  traced  for 


Figure  2.17  Geometry  for  the  determination  of  the  ray  separation, 

L. 
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exactly  the  same  length  of  time,  the  endpoints  of  these  rays  are  points  which  are 
in  phase  on  the  same  wavefront.  These  points,  {x\,zi)  and  {x2,Z2),  are  therefore 
connected  to  determine  the  ray  separation,  L. 


L  =  y{x2-  Xi)^  +  {Z2  -  Zxf 

The  ratio  of  to  as  in  Equation  2.35  may  finally  be  expressed  as 

P^  csin^oA^ 


(2.38) 


(2.3y) 


This  equation  is  only  practical  when  used  in  volume  cell  spreading  loss 
calculations.  For  closely  spaced  rays  ending  at  a  boundary,  either  the  surface 
or  the  bottom,  with  a  slight  difference  in  travel  time,  L  must  be  approximated 


L  =  Ar|  cos  (f>\ 


(2.40) 


where  0  is  the  angle  at  the  boundary  and  Ar  is  the  horizontal  separation  of  the 
two  rays  [2].  The  geometry  for  rays  ending  at  a  boundary  is  shown  in  Figure 


2.18.  Equation  2.40  may  be  substituted  into  2.39  to  obtain 

P^  _  csin<^oA<^ 

P^  Cof  Ar|  cos  <i>\ 


(2.41) 


which  is  valid  for  the  determination  of  spreading  loss  between  two  rays  at  a 
boundary. 

Since  acoustic  reciprocity  holds  true  for  a  non-isovelocity  environment  [7,16], 
the  two-way  transmission  loss  due  to  spreading  may  simply  be  expressed  as  twice 
the  one-way  spreading  loss  (in  dB)  [2].  This,  of  course,  assumes  that  the  path 
out  to  the  field  point  is  the  ssune  as  the  return  path.  If  the  path  to  the  scatterer 
were  different  from  the  return  path,  then  the  spreading  loss  associated  with  each 
path  would  have  to  be  determined  individually. 

Absorption  loss  in  a  horizontally  stratified  medium  may  not  be  assumed  to  be 
constant  for  all  raypaths  at  a  given  travel  time.  In  fact,  the  absorption  coefficient. 


Geometry  for  the  determination  of  the  ray  separation 
at  a  boundary. 
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a,  varies  with  temperature,  salinity,  and  pressure  throughout  the  medium.  In 
fresh  water,  the  absorption  coefficient  is  dependent  on  the  shear  viscosity,  /z,,  and 
the  volume  viscosity,  /it,,  of  the  water  as  well  as  on  the  square  of  the  frequency. 
The  fresh  water  absorption  coefficient,  ap-  may  be  expressed  as  follows  jl7j: 

8n^  /  3  \  ,2  ,  , 

In  sea  water,  however,  frequencies  below  about  lOOkHz  are  predominantly 
absorbed  through  an  ionic  relaxation  process.  This  moleculau-  dissociation- 
reassociation  process  occurs  with  fluctuations  of  pressure  2ts  an  acoustic  dis¬ 
turbance  passes.  A  finite  length  of  time,  the  relaxation  time,  is  required  for  this 
process  to  occur  [9j.  The  reciprocal  of  the  relaxation  time  is  referred  to  as  the 
relaxation  frequency  [9|.  Although  the  most  abundant  salt  found  in  sea  water  is 
sodium  chloride,  it  is  the  magnesium  sulfate,  MgSO^,  which  is  the  principle  cause 
of  absorption  in  the  range  from  about  IQkHz  to  lOOkHz  [18].  The  relaxation 
frequency  for  Mg  SO  a  may  be  expressed  as  [18] 

f(MgS04)  =  1-55  X  10^(r  +  273.1)  exp[-3052/(r  -I-  273.1)]Hz.  (2.43) 

Boric  acid,  B{OH)3,  dominates  the  absorption  below  about  IQkHz,  with  a 
relaxation  frequency  in  the  IkHz  region  [18]  which  may  be  expressed  as 

=  1-32  X  10^(r  -i-  273.1)  exp[-1700/(r  +  21Z.l)\Hz.  (2.44) 

Both  the  relaxation  processes  and  the  shear  and  volume  viscosities  must  be 
considered  in  order  to  obtain  the  absorption  coefficients  for  frequencies  up  to 
at  least  \QOkHz. 

Fisher  and  Simmons  [18]  developed  an  equation  for  the  absorption  coefficient 
in  sea  water  which  is  dependent  on  temperature,  pressure,  and  salinity.  The 
effects  of  the  relzocation  processes  of  MgSOA  and  B{OH)z  were  considered,  as 
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well  as  the  fresh  water  absorption  contribution.  The  salinity  was  assumed  to  be 
35  parts  per  thousand  (ppt)  and  the  pH  to  be  8.0.  The  results  of  this  work  are 
as  follows: 

a  =  8.686(Ai/i/V(/?  +  /")  +  A2P2/2/V(/2  +  f)  +  AaPa/'JdB/m.  (2.45) 
Where, 

Ai  =  (1.03  X  10-®  +  2.36  X  -  5.22  x  10-^^T^)sec/m, 

A2  =  (5.62  X  10-*  +  7.52  x  10“‘°r)sec/m, 

A3  =  (55.9  -  2.37T  +  4.77  x  -  3.48  x  x  IQ-^^scc^/m, 

P2  =  1  -  10.3  X  10-“*^  +  3.7  X  10-^P2, 
and, 

P3  =  1  -  3.84  X  lO'^P  +  7.57  X  10“®p2, 
for,  /  measured  in  Hz,  T  in  degrees  Celsius,  and  P  in  atmospheres.  The 
frequencies,  fi  and  /2,  axe  relaxation  frequencies  of  boric  acid  and  magnesium 
sulfate,  respectively. 

In  a  non-isovelocity  environment,  the  absorption  coefficient  varies  contin¬ 
uously  with  depth.  If  an  average  absorption  coefficient  is  assumed  constant 
throughout  each  constant  velocity  gradient  layer,  then  the  absorption  loss  along 
a  ray  in  any  layer  would  be  calculated  by  a  multiplication  of  the  length  of  the 
ray  segment  (arclength  of  a  circular  path)  by  the  average  absorption  coefficient 
(dB/m)  in  that  layer. 

ABSORPTION  LOSS^dB)  =  PATHLENGTH  x  a  (2.46) 

The  total  absorption  loss  at  the  end  of  a  ray  is  a  summation  of  the  losses  (in  dB) 
in  all  layers  through  which  the  ray  passes. 

When  acoustic  energy  reflects  from  a  boundziry,  such  as  the  ocean  surface 
or  bottom,  additional  transmission  losses  arise.  From  the  surface,  very  little 
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energy  is  2u:tually  transmitted  through  to  the  air  above,  as  a  pressure-release  type 
boundary  condition  exists.  The  impedance  mismatch  at  the  water-air  interface 
causes  the  sound  to  remain  contained  in  the  ocean.  In  addition  to  being  an 
excellent  reflector  of  sound,  the  ocean  surface  is  a  tremendous  scatterer  of  sound 
because  of  its  roughness  and  trapped  air  bubbles  or  biological  life  just  below 
it.  The  acoustic  energy  impinging  upon  the  surface  is  actually  scattered  in  all 
directions.  This  scattering  is  a  very  complex  function  of  wave  height,  wind 
speed,  rain,  grazing  angle,  and  frequency  [19,20].  An  additional  attenuation 
process  at  the  surface  is  acoustic  absorption  due  to  tiny  resonant  air  bubbles 
near  the  surface.  High  concentrations  of  trapped  bubbles  occur  near  the  surface 
due  to  environmental  conditions  (21).  A  relatively  high  amount  of  absorption 
results  from  small  (compared  to  a  wavelength)  bubbles  because  of  the  isothermal 
conditions  that  arise;  there  may  be  significant  heat  transfer  from  the  bubble  to 
the  surrounding  water.  Sound  propagation  is  usually  assumed  to  be  an  adiabatic 
process,  one  in  which  no  heat  transfer  occurs.  In  the  case  of  small  resonant 
bubbles,  however,  heat  will  transfer  from  the  air  to  the  water,  thus  yielding 
energy  absorption  [2].  Some  scattering  results  from  the  presence  of  bubbles,  but, 
the  absorbed  energy  is  much  greater  than  the  scattered  energy  [21|.The  sound 
energy  that  continues  forward  after  a  specular  reflection  is  lower  than  the  incident 
energy  primarily  as  a  result  of  the  energy  “lost”  to  scattering. 

Bottom  reflections  behave  quite  differently  than  surface  reflections.  In  many 
cases,  the  bottom  type  changes  gradually  with  depth.  The  impedance  also 
changes  gradually  from  that  of  the  ocean  water.  Thus,  some  of  the  impinging 
energy  is  free  to  be  transmitted  into  the  bottom  and  to  be  attenuated  due  to 
absorption.  The  incident  sound,  however,  may  additionally  be  scattered  in  all 
directions  as  well  as  reflected  forward.  In  other  cases,  the  bottom  may  be  very 
hard,  such  as  rock.  In  this  event,  little  energy  would  penetrate  and  be  absorbed 
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by  the  bottom.  Rather,  the  majority  of  incident  energy  would  reflect  forward. 
The  reflection,  absorption,  and  scattering  from  the  bottom  are  dependent  on  the 
bottom  type,  density,  speed  of  sound,  grazing  angle,  and  frequency.  Considerable 
work  has  been  done  on  the  formulation  of  expressions  for  bottom  scattering 
[22,23,24,25];  however,  there  is  still  a  substantial  amount  of  research  to  be  done 
in  this  field. 
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Chapter  3 

THE  INCORPORATION  OF  RAYTRACING  AND  MULTIPLE 
BOUNDARY  RETURNS  INTO  AN  ISOVELOCITY  MODEL 


3.1.  Isovelocity  Reverberation  Model 

A  computer  model  for  predicting  the  power  level  and  spectral  characteristics 
of  reverberation  has  been  developed  [1,26,27].  Effects  due  to  environmental 
conditions,  source  and  receiver  beam  patterns,  and  a  transmitted  waveform 
Eire  taken  into  account.  Although  reverberation  from  the  surface,  volume,  and 
bottom  are  considered  in  the  model,  no  forward  reflections  are  allowed.  Thus, 
whenever  any  portion  of  the  transmitted  energy  impinges  upon  either  the  surface 
or  the  bottom,  energy  is  assumed  to  be  backscattered  only;  i.e.,  forward-reflecting 
energy  is  neglected.  Forward-reflected  sound,  however,  may  be  very  significant. 
For  instance,  a  smooth  surface,  hard  bottom,  and  shallow  water  might  all  account 
for  significant  reverberant  returns  after  reflections  occur. 

Also  assumed  in  this  model  is  an  isovelocity  profile;  i.e.,  the  speed  of  sound 
has  been  assumed  to  be  constant  everywhere  throughout  the  medium.  The 
geometry  of  the  problem  has  been  greatly  simplified  through  this  assumption. 
Refraction  of  sound  is  non-existent  under  isovelocity  conditions,  and  the  acoustic 
energy  spreads  out  spherically  from  the  source.  The  isovelocity  approximation, 
however,  may  not  be  a  good  one,  especially  over  long  distances  or  around  regions 
in  which  strong  thermal  or  salinity  gradients  exist. 

3.1.1.  Introduction  to  REVMOD 

The  base  isovelocity  reverberation  model  to  be  modified,  REVMOD  (REVer- 
beration  MODel)  [1,26,27],  is  a  package  consisting  of  five  FORTRAN  computer 
programs.  In  this  model,  the  ocean  is  divided  into  a  finite  number  of  cells, 
Nf  at  a  distance  from  the  source  which  is  half  of  the  round-trip  distance  from 
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the  source  back  to  the  receiver.  These  cells  are  the  scatterers  whose  reverberant 
returns  are  summed  to  obtain  a  composite  reverberant  spectrum  or  a  total  power 
level.  For  each  distance  of  interest,  a  scattering  function  is  accumulated.  This 
scattering  function  is  then  convolved  with  the  power  spectrum  of  the  transmitted 
waveform  [1,4]  to  obtain  the  total  reverberation  power  spectrum  at  any  distance. 
See  Figure  3.1  for  a  flow  diagram  of  how  REVMOD  computes  the  reverberation 
power  level  and  spectrum.  The  purpose  of  this  study  is  to  incorporate  raytracing 
and  multiple  boundary  returns  into  the  model.  The  section  which  produces 
the  scattering  function,  SCATTER,  is  of  primary  interest.  Following  is  a  brief 
description  of  each  of  the  progranos  that  make  up  the  unmodified  version  of 
REVMOD,  a  synopsis  of  the  major  reverberation  modeling  program,  SCATTER, 
a  discussion  of  the  capabilities  of  REVMOD  and  some  sample  output. 

In  the  first  program  section  of  the  REVMOD  series,  BEAMS,  files  containing 
desired  transmitting  and  receiving  beam  patterns  for  later  use  in  the  reverber¬ 
ation  model  are  created.  Both  beams  are  assumed  to  be  circularly  symmetric. 
Beam  pattern  weightings,  normalized  to  one  at  zero  degrees,  are  input  from  zero 
to  ninety  degrees.  Each  beam  may  be  pointed  in  an  arbitrary  direction  through 
the  specification  of  a  beam  vector,  (0,^),  where  ^  is  an  azimuthal  angle  and  <!>  is 
an  elevation  mgle.  An  option  has  also  been  made  available  for  the  utilization  of 
two-dimensional  beam  patterns,  but  the  option  is  not  exercised  in  this  study. 

All  of  the  information  needed  for  the  scattering  function  calculations  is 
accepted  in  the  input  program  of  REVMOD,  INPUT.  Included  as  input  data  are 
the  source  speed  and  depth,  the  environmental  conditions  of  the  ocean,  the  pulse 
length  and  frequency,  the  acoustic  beam  patterns,  and  the  grid  cell  dimensions. 
See  Appendix  C  for  a  complete  list  of  the  INPUT  variables. 

The  scattering  function  of  the  ocean  for  a  specified  distance,  or  series  of 
distances  is  computed  in  the  third  program,  SCATTER.  In  this  program,  a  grid  is 


Figure  3.1  REVMOD  flow  diagram. 
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used  to  divide  the  ensonified  boundaries  and  volume  of  the  ocean  into  cells.  Each 
cell  is  individually  processed  for  its  contributions  to  the  backscattered  spectrum 
or  power  gain.  Taken  into  account  are  spectral  spreading  and  shifting,  and 
attenuation  factors  resulting  from  the  environmental  conditions,  the  geometry, 
and  the  beam  patterns.  The  SONAR  (SOund  NAvigation  Ranging)  equation 
is  solved  for  each  cell  to  determine  the  relative  backscattering  level  [26] .  The 
actual  transmitted  power  level  is  not  needed  at  this  time  as  it  is  assumed  that  the 
scattering  function  is  independent  of  the  source  level.  In  addition  to  calculating 
the  total  spectrum,  SCATTER  has  the  capability  of  separating  the  surface, 
volume,  and  bottom  contributions  so  that  each  may  be  studied  individually. 

Until  this  point  in  the  reverberation  model,  the  source  level  and  signal 
envelope  have  not  been  of  concern.  The  scattering  function  which  characterizes 
the  ocean  environment  is  independent  of  these  parameters.  In  the  fourth  prograun 
of  the  REVMOD  series,  TRANSMIT,  the  normalized  energy  spectrum  of  the 
transmit  pulse  is  convolved  with  the  scattering  function  produced  by  SCATTER 
to  obtain  the  reverberation  power  spectrum  at  the  receiver  [26]. 

3.1.2.  Description  of  SCATTER 

The  incorporation  of  raytracing  and  multiple  boundary 
returns  into  REVMOD  involved  a  rewriting  of  program  SCATTER.  Although 
some  SCATTER  subroutines  have  been  kept  intact,  many  have  been  altered 
or  excluded  altogether.  In  this  section,  the  functions  of  all  pertinent  SCATTER 
subroutines,  as  they  existed  before  any  changes  had  been  made,  are  described.  It 
will  serve  as  a  reference  later  when  the  raytracing  and  multiple  boundary  return 
algorithms  are  discussed.  A  flow  diagram  for  SCATTER  is  shown  in  Figure  3.2a. 
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Figure  3.2  Flow  diagram  for  SCATTER. 
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b.  Spectral  density  section  A. 
Figure  3.2  Continued. 
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The  section  of  this  flow  diagram  which  determines  the  spectral  density  is  shown 
in  Figure  3.2b.  All  of  the  subroutines  which  appear  in  Figures  3.2a  and  3.26  are 
discussed  either  in  this  section  or  in  Appendix  D. 

Provided  in  Appendix  D.l  is  a  brief  description  of  the  initialization  subrou¬ 
tines  found  in  SCATTER.  Once  all  input  parameters  are  initialized,  the  SONAR 
equation  is  solved  in  three  steps.  First,  subroutine  SEQNl  is  called  once  for 
each  distance,  at  which  time  a  partial  SONAR  equation  solution  is  calculated. 
Determined  in  SEQNl,  is  the  two-way  transmission  loss,  which,  in  this  isovelocity 
model,  is  the  same  for  all  cells  at  one  distance,  regardless  of  the  scattering 
mechanism,  surface,  volume,  or  bottom.  Only  absorption  and  spreading  losses 
are  considered  in  this  section. 

Secondly,  in  subroutine  SEQN2,  a  partial  SONAR  equation  for  each  type 
of  reverberation,  only  one  type  per  cell,  is  calculated.  See  Appendix  D.2  for  a 
brief  outline  of  how  the  different  reverberation  types  are  distinguished.  SEQN2  is 
called  once  for  each  row  of  cells  as  the  backscattering  mechmism,  surface,  volume, 
or  bottom,  is  constant  throughout  any  one  row.  In  SEQN2,  a  backscattering 
coefficient  is  multiplied  by  the  ensonified  area  of  the  scatterer  to  obtain  the 
relative  backscattered  intensity  level.  If  a  volume  cell  is  being  processed,  the 
backscattering  coefficient  is  multiplied  by  the  volume  of  the  scattering  element. 
The  partial  SONAR  equation  computed  in  SEQN2  is  multiplied  by  the  value 
computed  in  SEQNl  to  obtain  a  new  partial  SONAR  equation  value. 

Finally,  the  SONAR  equation  is  completed  for  each  cell  in  subroutine  SEQN3. 
Here,  the  partial  SONAR  equation  from  SEQN2,  which  accounts  for  spreading, 
absorption,  and  backscattering  level,  is  multiplied  by  the  composite  beam  pattern 
attenuation.  Both  transmitting  and  receiving  beam  patterns  are  taken  into 
account.  See  Appendix  D.3  for  a  description  of  subroutine  BPATRM,  in  which 
the  beam  attenuations  are  interpolated.  Once  SEQN3  has  been  called  for  a  cell. 
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that  cell’s  contribution  to  the  power  gain  is  complete.  However,  if  a  reverberation 
spectrum  is  desired,  some  additional  operations  must  be  performed  on  that  cell. 
See  Appendix  DA  for  the  descriptions  of  subroutines  ZPDF,  DOPP,  CNVOLl, 
SWPDF,  WTSCAT,  and  RINIT2  and  RINIT3.  These  include  operations  which 
determine  the  Doppler  shift  due  to  source  motion,  current  layers  and  wind  speed, 
and  the  spectral  spreading  due  to  random  scatterer  motion. 

3.1.3.  Capabilities  and  Output  from  SCATTER 

Program  SCATTER  has  been  designed  to  produce  two  types  of  output.  The 
first  type  is  a  plot  of  the  power  gain  of  the  scattering  function  versus  time.  The 
variation  in  the  received  reverberation  level  with  time  may  be  examined  through 
this  feature.  Plots  of  this  type  are  illustrated  in  Figure  3.3.  The  individual 
contributions  to  the  scattering  function  from  the  surface,  voliune,  and  bottom 
axe  shown  in  addition  to  the  combined  total.  The  scenario  corresponding  to 
the  plot  in  Figure  3.3  is  shown  in  Figure  3.4.  The  input  data  used  may  be 
found  in  Appendix  E.l.  The  first  sxuface  return  is  marked  by  point  S  in  Figures 
3.3  and  3.4.  It  is  seen  to  appear  at  a  time  of  0.40  seconds  after  transmitting. 
The  first  bottom  return  is,  likewise,  maurked  by  point  B  in  Figures  3.3  and  3.4. 
This  return  is  seen  first  at  a  time  of  0.27  seconds  after  transmitting.  Factors 
such  as  the  scattering  strengths,  absorption  coefficient,  and  beam  patterns  are 
not  necessary  to  be  known  in  the  understanding  of  thb  introductory  example. 
However,  a  major  role  is  played  by  these  factors  in  the  shaping  of  the  general 
curves  of  Figure  3.3. 

The  second  type  of  output  that  SCATTER  has  been  designed  to  yield  is 
the  spectral  density  of  the  scattering  function  for  one  time.  The  spreading  and 
shifting  in  frequency  of  the  transmitted  signal  is  dependent  on  both  source  and 
scatterer  motion.  The  source  motion  has  been  assumed  to  be  a  constant  velocity 
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Figure  3.3  Output  from  SCATTER.  Power  gain  versus  time  plots 
for  the  surface,  volume,  bottom  and  combined  total. 
Ocean  depth:  500  m. 
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in  the  horizontal  plane.  The  scatterer  motion  is  modeled  as  two  components, 
a  zero-mean  Gaussian  velocity  distribution  and  an  average  velocity  vector. 
The  former  is  the  cause  for  a  spreading  of  frequency  about  some  transmitted 
frequency.  This  is  actually  modeled  as  a  random  Doppler  shifting.  The  latter, 
along  with  the  source  motion,  is  the  cause  of  the  mean  Doppler  shifting  of  the 
transmitted  signal.  Some  basic  scattering  functions  are  shown  in  Figure  3.5. 
The  velocity  for  the  source  was  chosen  to  be  5  m/s  and  the  elapsed  time  from 
transmitting  to  receiving  was  set  as  0.48  seconds.  Aside  from  the  fact  that  the 
source  has  been  assigned  a  velocity,  the  scenario  used  for  Figure  3.5  is  identical  to 
that  shown  in  Figure  3.4.  The  inputs  aire  shown  in  Appendix  E.2.  The  acoustic 
beams  used  in  these  examples  are  omni-directional  in  the  vertical,  direction 
and  ten  degrees  wide  in  the  horizontal,  d,  direction  with  no  sidelobes  present. 
These  beaims  have  been  chosen  to  place  an  emphasis  on  the  effects  from  the 
surface  and  bottom  returns,  and  thus,  to  best  illustrate  the  differences  of  the 
modified  scattering  function  from  those  of  the  original  model. 

Before  a  disc\ission  of  the  significance  of  Figure  3.5,  a  brief  explanation  of 
the  Doppler-shifting  mechanisms  is  necessauy.  In  SCATTER,  a  frequency  offset 
has  been  introduced,  such  that  a  stationary  scatterer,  which  lies  in  line  with 
the  vector  of  source  motion,  produces  zero  Doppler  shift.  This  frequency  offset 
has  an  effect  on  each  cell.  The  cells  which  lie  ninety  degrees  off  of  the  source 
motion  axis  account  for  the  greatest  amount  of  Doppler  shift  due  to  the  offset, 
although  their  motion  with  respect  to  the  source  is  minimal.  This  shift  occurs  in 
the  negative  Doppler  region  as  a  result  of  the  frequency  offset.  When  the  source 
is  stationary,  no  frequency  offset  is  introduced  and  the  Doppler  from  each  cell 
is  related  only  to  the  motion  of  the  cell  with  respect  to  the  stationary  source. 
In  Figure  3.5,  both  a  spreading  and  shifting  of  frequency  can  be  observed.  The 
frequency  spreading  that  is  present  in  ail  ceils  is  due  to  the  standard  deviation 
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of  the  scatterer  speed  distributions.  The  shifting  here  is  a  result  of  the  source 
motion  only.  All  scattering  elements  have  a  mean  velocity  of  zero. 

Figure  3.5  shows  a  single  surface  return,  a  single  bottom  return  and  a  volume 
reverberation  return.  The  surface  return  is  the  more  negatively  Doppler  shifted 
one,  as  the  grzizing  angle  at  the  surface  is  greater  than  the  grazing  angle  at 
the  bottom.  The  volume  return  in  Figure  3.5  mainly  underlies  the  surfa.ce  and 
the  bottom  contributions.  However,  some  positive  frequency  spreading  from  the 
volume  may  be  observed.  This  is  due  to  the  random  motion  of  those  volume 
elements  which  lie  on,  or  very  close  to,  the  axis  of  source  motion. 

3.2.  Modified  Reverberation  Model 

A  modified  version  of  SCATTER  has  been  written  to  construct  a  scattering 
function  for  a  non-isovelocity  environment  in  which  reflections  may  occur  at 
the  boundaries.  This  new  program,  called  RAYSCT  (RAYtracing  version 
of  SCaTter),  has  been  designed  to  use  the  same  approach  as  SCATTER  in 
producing  a  scattering  function.  The  concept  of  a  grid  of  cells  is  maintained, 
as  well  as  the  program’s  ability  to  separate  the  surface,  volume,  and  bottom 
contributions  to  reverberation  from  one  another.  The  scattering  function,  which 
is  produced  in  RAYSCT,  is  compatible  with  the  TRANSMIT  section  of  the 
REVMOD  series.  In  other  words,  the  output  scattering  function  from  RAYSCT 
may  be  convolved  with  the  power  spectrum  of  any  of  the  transmit  signal  types 
found  in  TRANSMIT  to  obtain  the  total  reverberation  power  spectrum. 

3.2.1.  Considerations  for  the  Non-isovelocity  Model 

Although  the  general  format  for  the  determination  of  the  scattering  function 
has  remained  unchanged,  there  have  been  many  aspects  of  the  problem  which  had 
to  have  been  considered  carefully.  For  one,  the  grid  cells  in  the  modified  model 
are  not  consistent  in  size.  They  are  also  summed  over  a  contour  of  constant  time 
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as  opposed  to  the  spherical  contour  of  constant  distance  found  in  SCATTER. 
Another  complication  is  that  the  cells,  which  fall  on  the  surface  and  bottom 
boundaries  at  a  certain  time,  are  not  as  clearly  defined  due  to  both  the  non- 
isovelocity  environment  and  the  allowance  of  reflections.  In  fact,  there  may  be 
a  number  of  different  patches  of  the  surface  and  bottom  which  are  ensonified 
simultaneously.  Also,  as  opposed  to  the  isovelocity  model,  SCATTER,  every 
single  cell  is  a  volume  reverberation  cell,  which  may  or  may  not  ensonify  a 
boundary  as  well.  Additionally,  there  are  some  complicated  Doppler  situations 
which  may  arise  due  to  varying  speeds  of  sound  throughout  the  medium. 
Furthermore,  the  backscattered  Gaussian  frequency  distribution  is  varied  from 
cell  to  cell  due  to  the  differences  in  sound  speed.  Finally,  as  reflections  were 
not  permitted  in  SCATTER,  a  procedure  for  subtracting  energy  which  is  lost 
to  interactions  with  the  surface  and  bottom  boundaries  had  to  be  developed. 
Each  of  these  factors  plays  an  important  role  in  the  development  of  the  modified 
scattering  function  program,  RAYSCT. 

The  grid  of  cells  in  RAYSCT  was  not  able  to  be  defined  as  uniformly  as 
it  had  been  in  SCATTER.  The  critical  angles  that  were  defined  in  Chapter  2 
are  used  to  define  different  sectors  of  which  the  entire  grid  is  made.  Figure  3.6 
illustrates  a  set  of  critical  angles  being  used  to  divide  the  ensonified  medium 
into  what  may  be  defined  as  critical  sectors.  In  RAYSCT,  these  sectors  are  each 
subdivided  into  equal  vertical  angular  increments,  such  that  the  sector  contains 
the  smallest  whole  number  of  increments  which  are,  in  size,  less  than  or  equal 
to  a  predetermined  maximum  grid  cell  height.  For  example,  if  a  critical  sector 
were  to  have  a  vertical  spread  of  11®,  and  the  maximum  grid  cell  height  allowed 
by  a  user  was  2°,  then  the  sector  would  be  subdivided  into  6  cells,  each  having  a 
vertical  spread  of  1.83®.  Five  subdivisions  would  have  caused  the  cells  to  be  too 
large  and  seven  would  have  been  more  than  were  necessary.  Since  the  critical 


sectors  may  vary  in  size,  the  cells  in  one  critical  sector  would  not  necessarily  be 
the  same  size  as  the  cells  in  another  critical  sector.  No  additional  problems  are 
brought  about  in  the  event  of  differently  sized  cells  because  the  backscattered 
level  from  each  cell  is  dependent  on  the  angular  separation  of  the  rays  to  that 
cell. 

In  program  SCATTER,  the  ensonified  medium  is  divided  into  three  regions; 
surface,  volume,  and  bottom.  In  general,  the  cells  with  the  greatest  elevation  are 
considered  to  yield  surface  reverberation.  The  middle  region  is  comprised  of  cells 
which  yield  volume  reverberation  and  cells  with  the  least  elevation  yield  bottom 
reverberation.  Figure  3.7  is  a  general  case  of  a  grid  of  ceils.  All  three  types 
of  reverberation,  however,  might  not  exist  simultaneously  at  certain  times.  In 
Figure  3.7,  the  boundaries  dividing  each  type  of  cell  are  well  defined.  In  a  non> 
isovelocity  environment,  the  regions  of  boundary  ensonification  are  not  as  well 
defined.  Once  multiple  boundary  returns  are  allowed  via  reflections,  there  might 
be  any  number  of  groups  of  cells  which  ensonify  the  surface  or  the  bottom  at 
the  same  instant.  Figure  2.3  is  an  example  of  different  regions  of  each  boundary 
being  ensonified  simultaneously.  Also,  a  reverberant  return  from  a  boundary 
at  one  time  may  vanish  completely  at  some  later  time  due  to  effects  from  the 
sound  velocity  profile.  Figure  3.8  illustrates  how  the  acoustic  energy  is  refracted 
downward  and  away  from  the  surface  by  a  negative  velocity  gradient,  such  that, 
at  the  longer  time,  the  direct  ray  path  from  the  source  to  the  surface  has  vanished. 

Another  matter  which  was  considered  in  the  modification  of  the  scattering 
function  computing  program  is  the  fact  that  all  boundary  cells,  those  falling  on 
either  the  surface  or  the  ’  ottom,  will  yield  two  types  of  reverberation.  One  type 
is  from  the  particular  boundary  of  concern  and  the  other  is  volume  reverberation. 
In  Figure  3.9,  all  surf2u:e  cells  are  accompanied  by  a  fraction  of  a  volume  cell. 
The  additional  volume  cell  fractions,  which  are  crossed  out  in  Figure  3.9,  are  not 


71 


Figure  3.7  Differentiation  between  surface,  volume,  and  bottom 
cells.  Isovelocity  environment,  reflections  are  disal¬ 
lowed. 


Figure  3.8 


A  direct  surface  return  vanishes  at  some  later  time  in 
the  presence  of  a  negative  velocity  gradient. 
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Figure  3.9  Celia  marked  with  an  “X"  have  been  assumed  to  yield 
no  volume  reverberation  in  SCATTER. 
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taken  into  account  in  SCATTER.  Each  cell  has  been  assumed  to  yield  one  type 
of  reverberation  only.  However,  in  program  RAYSCT,  the  fractions  of  volume 
cells  which  fall  both  before  and  after  a  boundary  return  have  been  considered. 
Figure  3.10  illustrates  a  surface  cell  that  is  located  inside  a  volume  cell.  Each 
part  of  the  volume  cell  must  be  weighted  differently,  not  only  because  of  the 
differences  in  size,  but  because  the  portion  after  the  reflection  would  have  lost 
some  additional  energy  due  to  the  inter2u:tion  with  the  surface. 

In  order  for  a  cell  to  be  divided  at  a  boundary,  at  least  two  ray  segment 
lengths  must  be  known.  First,  the  total  grid  cell  depth  in  meters  should  be 
found.  This  is  determined  by  subtracting  the  ray  to  the  center  elevation  of  the 
front  of  the  cell  from  the  ray  to  the  center  elevation  of  the  back  of  the  cell. 
Next,  the  length  of  either  portion  of  the  divided  cell  must  be  determined.  For 
thb  dimension  to  be  found,  however,  the  length  of  the  ray  from  the  source  to 
the  boundary  must  first  be  determined.  For  use  in  RAYSCT,  a  set  of  ray  state 
numbers  similar  to  those  suggested  by  Roeckel  [14)  in  section  2.2.5  have  been 
developed.  KA  has  been  designated  as  the  number  of  times  a  ray  passes  through 
the  region  above  the  source  and  KB  has  been  designated  as  the  number  of  times 
a  ray  passes  through  the  region  below  the  source. 

A  raytracing  scheme  has  been  employed  such  that,  if  an  initial  angle, 
were  given,  the  time,  for  the  ray  to  travel  from  the  source  to  the  surface,  or 
until  a  turning  point  above  the  source  is  reached,  could  be  calculated.  If  this  ray 
were  to  continue,  it  would  reflect  from  the  surface,  or  turn,  and  eventually  pass 
through  the  depth  of  the  source  at  an  angle  of  —  due  to  symmetry.  Thus,  the 
raytr2M:ing  scheme  would  also  calculate  the  time,  of  the  ray,  transmitted  at 
to  travel  from  the  source  to  the  bottom,  or  until  it  reaches  a  turning  point 
below  the  source.  If  these  partial  travel  times,  as  well  as  the  state  numbers,  KA 
and  KB,  are  known,  then  the  time  for  the  ray  to  arrive  at  the  boundary  may  be 
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determined. 

^BOUNDARY  =  X  X  (3-1) 

Figure  3.11  shows  a  ray  from  the  source  to  the  surface  after  one  surface  and 
one  bottom  reflection.  For  this  case,  KA  is  equal  to  3  and  KB  is  equal  to  2.  If 
the  travel  times  and  t g  were  known,  the  time  for  the  ray  to  come  into  contact 
with  the  surface  the  second  time  would  be 


isURFACB  =  3  X  +  2  X  . 

Once  the  time  for  the  ray  to  travel  to  the  boundary,  tgou NDARY  i  known,  the 
distance  along  the  ray  from  the  source  to  the  boundary  is  determined  by  actually 
tracing  the  ray  for  this  amount  of  time  and  accumulating  the  ray  segment  lengths 
through  each  layer. 

Figure  3.12  shows  the  case  of  a  volume  reverberation  cell  which  is  split  by 
a  surface  reflection.  In  Figure  3.12,  the  total  cell  length  is  represented  by  ray 
segment  AC,  where 

AC  =  ^  -JA.  (3.2) 

Since  SB  is  known  through  the  t  BOUND  ARY  calculation, 

TB  =  'SB  -  JA.  (3.3) 


Similarly, 


BC  =  SC  -  SB. 


(3.4) 


With  ail  ray  segments  known,  the  weighting  factor,  Wi,  for  the  portion  of  the 
cell  before  the  reflection  in  Figure  3.12  would  be 


AB 

w  ■ 


(3.5) 
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SURFACE 


I 


Figure  3.11 


Illustration  of  a  ray  ending  at  the  surface  of  the  oceami. 
The  ray  state  numbers  are,  KA  =  3  and  KB  =  2. 
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BOTTOM 


Figure  3.12 


Geometry  for  the  division  of  a  volume  cell  at  the 
surface. 


It  follows  that  the  weighting  factor,  IV^,  for  the  portion  of  the  cell  after  the 
reflection  would  be 


BC 

Wi  = 

AC 


(3.6) 


such  that 


Wi+W2  =  1.0. 


(3.7) 


Thus,  energy  is  neither  added  to  nor  subtracted  from  the  cell  as  it  is  divided  and 
its  parts  processed  individually. 

There  is  a  problem  which  may  arise  in  the  determination  of  whether  a  volume 
cell  is  split  by  a  boundary  reflection.  A  part  of  the  cell  may  have  undergone  a 
reflection,  but  the  center  ray  of  that  cell  might  not  yet  have  reached  the  boundary. 
Figure  3.13  illustrates  this  case.  If  this  volume  cell  were  to  be  processed  as  was 
just  outlined,  the  entire  cell  would  be  weighted  uniformly,  as  the  center  ray  has 
not  reflected  at  the  boundary.  To  avoid  this  potential  trouble  spot,  an  additional 
set  of  angles,  which  will  be  referred  to  as  surface  and  bottom-solution  sector 
angles,  have  been  added  to  the  list  of  critical  angles.  Included  in  this  set  are  all 
angles  ending  at  either  boundary  at  a  time  corresponding  to  either  the  front  or 
the  back  of  the  scattering  volume.  In  Figure  3.14,  a  scattering  volume  is  shown 
to  contact  the  surface  by  a  direct  path  from  the  source.  A  surface-solution  sector 
is  defined  by  the  two  rays  shown.  This  solution  sector  may  be  made  up  of  many 
grid  cells.  Surface  reverberation  is  represented  by  all  cells  falling  between  the  two 
rays  in  Figure  3.14.  Each  of  these  cells,  though,  is  also  a  volume  reverberation 
cell  which  is  divided  at  the  surface.  No  cell  which  lies  outside  the  surface-solution 
sector  touches  the  surface  at  any  time.  With  the  solution  sector  defined  as  such, 
the  center  ray  of  any  volume  cell  which  falls  entirely  within  the  surface-solution 
sector,  is  reflected  from  the  surface,  thus,  making  it  impossible  for  the  situation 
shown  in  Figure  3.13  to  exist.  A  volume  cell  which  has  been  split  by  the  surface  is 


Figure  3.13  Partial  reflection  of  a  volume  ceil.  In  this  case,  the 
center  ray  of  the  cell,  AC^  may  not  be  used  to  indicate 
whether  a  reflection  has  occurred. 


Figure  3.14  The  shaded  area  is  a  surface  solution  sector. 
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illustrated  in  Figure  3.15.  Examples  of  volume  cells  which  fall  outside  the  surface- 
solution  sector  aire  shown  in  Figures  3.16  and  3.17.  An  additional  surface-solution 
sector  is  produced  each  time  a  surface  reflection  occurs.  Existing  in  the  same 
manner  as  surface-solution  sectors,  of  course,  are  bottom-solution  sectors. 

In  RAYSCT,  a  complete  set  of  boundary-solution  angles,  resulting  from  all 
surface  and  bottom  reflections,  is  accumulated  and  added  to  the  list  of  critical 
angles.  However,  after  the  reverberation  at  one  time  has  been  computed,  and  the 
next  time  increment  is  considered,  all  of  the  boundary-solution  angles  change. 
The  angles  intersecting  the  surface  and  bottom  are  generally  more  shallow  for 
a  longer  time  interval  between  transmitting  and  receiving.  Thus,  a  new  set  of 
boundary-solution  angles  must  be  determined  for  each  time  increment.  The  list 
of  critical  angles,  to  which  these  axe  added,  however,  remains  constant.  Their 
only  dependencies  axe  on  the  source  location  and  the  sound  velocity  profile. 

When  a  non-isovelocity  environment  is  present,  the  Doppler  shift  component 
for  a  moving  scatterer  is  not  determined  by  using  the  same  parameters  that  au^e 
used  for  a  moving  source.  The  speed  of  sound,  as  well  as  the  angle  of  the  ray,  may 
vaxy  from  the  source  location  to  the  scatterer  location.  It  turns  out,  however, 
that,  due  to  Snell’s  law,  the  effect  of  a  different  angle  at  the  scatterer  is  negated 
by  the  effect  of  a  different  speed  of  sound,  with  a  reference  to  the  transmitted 
angle  and  the  speed  of  sound  at  the  source.  Appendix  F  offers  a  more  detailed 
explanation.  Hence,  in  the  case  of  a  monostatic  return,  the  speed  of  sound  and 
trauismitted  angle  at  the  source  may  actually  be  used  along  with  the  velocity  of  a 
scatterer  to  determine  the  Doppler  shift  caused  by  that  scatterer.  In  RAYSCT, 
the  speed  of  sound  and  incident  angle  at  the  scatterer  have  been  used,  first,  so 
that  the  option  of  implementing  bistatic  scatter  to  the  problem  may  be  left  open 
for  future  research,  and  secondly,  simply  because  the  parameters  needed  were 
readily  accessible  from  the  raytracing  programs. 


Figure  3.17  The  volume  cell  shown  falls  below  the  surface-solution 
sector. 
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A  similar  case  has  arisen  concerning  the  backscattered  Gaussian  frequency 
distribution  for  each  ceil.  Although  all  ceils  of  the  same  type,  surface,  volume, 
or  bottom,  have  been  assumed  to  have  the  same  random  velocity  distribution, 
the  Gaussian  frequency  distributions,  which  are  a  direct  result  of  the  velocity 
distribution,  may  vary  from  cell  to  cell.  This  is  due  to  the  varying  speed  of 
sound  throughout  the  medium.  The  Gaussian  frequency  distribution  is  calculated 
as  a  random  Doppler  shifting  through  the  use  of  the  Gaussian  scatterer  speed 
distribution,  the  transmitted  frequency,  and  the  speed  of  sound  at  the  cell. 
Since,  in  RAYSCT,  the  speed  of  sound  at  the  scatterer  may  vary,  so  may  the 
Gaussian  frequency  distribution.  This  distribution  is  calculated  once  for  each 
horizontal  row  of  cells  in  RAYSCT.  Actually,  because  the  speed  of  sound  is 
constant  everywhere  along  the  surface  amd  the  bottom,  only  volume  cells  can 
have  a  varying  Gaussian  frequency  distribution. 

Since  reflections  are  not  permitted  in  the  isovelocity  reverberation  model, 
some  assumptions  had  to  be  made  in  its  modification.  First,  specular  reflections 
were  assumed,  in  which  the  angle  of  incidence  is  equal  to  the  angle  of  reflection. 
It  was  assumed  that  no  frequency  spreading  occurs  during  a  reflection;  the 
only  frequency  spreading  allowed  is  that  which  occurs  during  the  backscattering 
process.  Another  assumption  that  was  made  is  that  losses  due  to  surface  and 
bottom  reflections  are  independent  of  both  frequency  and  grazing  angle.  The 
incorporation  of  complicated  reflection  loss  functions  into  RAYSCT  in  the  future 
is  anticipated.  The  code  has  been  written  for  the  program  to  conveniently  accept 
these  functions.  At  the  present,  though,  two  loss  factors,  one  for  the  surface  and 
one  for  the  bottom,  are  used  as  input  variables  to  the  modified  program. 
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Figure  3.18  Program  RAYSCT. 
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3.2.2.  Modified  Algorithm  for  the  Scattering  Function  Calculation 

Several  raytracing  subroutines  have  been  written  for  the  modified  version  of 
the  scattering  function  program,  RAYSCT.  Figure  3.18  is  a  flowchart  of  the  basic 
operation  of  this  program.  A  hierarchy  diagram  of  the  subroutines  called  from 
RAYSCT  to  produce  the  modified  scattering  function  is  shown  in  Appendix 
G.  Some  of  these  routines  had  been  found  originally  in  SCATTER,  some  are 
based  on  subroutines  from  SCATTER,  and  some  are  new  raytracing  subroutines 
written  primarily  for  the  modification  of  SCATTER.  The  asterisks  in  Appendix 
G  signify  subroutines  from  SCATTER  which  were  left  unchanged  in  the  modified 
program.  In  this  section,  the  functions  of  the  new  subroutines  which  appear  in 
RAYSCT  are  discxissed. 

The  input  to  the  new  program,  RAYSCT,  is  quite  similar  to  the  input  for 
the  original  program,  SCATTER.  However,  there  are  some  important  differences. 
Program  RAYSCT  accepts  a  sound  velocity  profile  (SVP),  instead  of  a  constant 
speed  of  sound.  This  SVP  consists  of  up  to  fifty  depth-velocity  point  pairs.  The 
ocean  depth  need  not  be  input  since  it  is  predetermined  by  the  maocimum  depth 
of  the  SVP.  An  absorption  coefficient  is  calculated  for  each  constant  velocity 
gradient  layer  in  the  SVP.  Power  loss  factors  for  surface  and  bottom  refiections 
are  added  to  the  list  of  input  variables.  The  total  spread  in  elevation  of  the 
grid  is  always  assumed  to  be  one  hundred  eighty  degrees.  An  upper  limit  on  the 
vertical  cell  size  is  an  input  variable  as  the  vertical  cell  size  may  vary.  Finally, 
the  input  variables  for  distances  in  meters  have  been  changed  to  times  in  seconds 
where  appropriate. 

Several  of  the  raytracing  routines,  which  had  been  written  for  the  modified 
program,  initialize  certain  parameters.  These  parameters  are  used  to  set  up  a 
new  grid  for  the  incremental  scattering  function  accumulation.  In  the  subroutine 
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CRITICAL-ANGLES  (Appendix  a  complete  set  of  critical  angles  between 
plus  ai'  1  minus  ninety  degrees,  with  a  reference  of  zero  degrees  as  horizontal,  are 
calculated.  These  are  used  to  define  the  critical  sectors,  regions  throughout 
which  ail  rays  behave  similarly.  The  critical  sectors  are  actually  formed  in 
RAYSET  (Appendix  H.2)  from  the  list  of  critical  angles.  In  the  subroutine 
RAYBOUND  (Appendix  H.3),  these  critical  sectors  are  searched  for  all  sets  of 
rays  which  hit  a  surface  or  bottom  boundary  within  the  scattering  volume.  Each 
set  of  rays  determined  in  RAYBOUND  is  designated  as  either  a  surface  or  a 
bottom-solution  sector.  Subroutines  TIMEBOUND,  BROKEN,  ITERATION, 
and  SOLBND  (Appendices  HA,  H.5,  H.6,  and  H.7  respectively)  are  each  called 
to  iterate  upon  a  solution  for  the  angle  of  the  ray  which  reaches  a  particular 
boundary  at  a  specified  time.  In  the  last  subroutine  of  the  set  up  procedure, 
VOLANGS  (Appendix  H.8),  a  list  of  volume  angles  is  compiled  by  the  addition 
of  all  surface  and  bottom-solution  angles  to  the  set  of  critical  angles.  These 
angles  are  called  volume  angles  beca;ise  they  are  used  to  divide  the  ensonified 
volume  into  sectors.  The  sectors  are  later  subdivided  into  vertical  grid  cells. 

Once  all  of  the  sectors  have  been  defined  for  a  particular  time,  the  scattering 
function  accumulations  are  begun.  Whereas  one  loop  to  sweep  through  all  cells 
from  the  surface  to  the  bottom  is  used  in  SCATTER,  three  different  subroutines 
arc  called  to  perform  this  function  in  RAYSCT.  These  subroutines  are  called 
SURSEQ  (SURface  Sonar  EQuation),  VOLSSQ  (VOLume  Sonar  EQuation), 
and  BOTSEQ  (BOTtom  Sonar  EQuation).  Flowcharts  for  these  subroutines 
are  located  in  Appendices  H.9,  H.ll,  and  H.12  respectively.  Appendix  H.IO 
is  a  flowchart  for  SPECTRAL  DENSITY  SECTION  B,  which  is  referred  to  in 
Appendix  H.9. 

In  the  subroutine  SURSEQ,  every  surface-solution  sector  is  processed 
individually.  Each  sector  is  divided  into  equal-sized  vertical  angular  increments 
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whic}-  are  smaller  than  or  equal  to  an  input  maximum  vertical  cell  size.  The 
horizontal  angular  cell  size,  Ad,  is  constant  for  all  cells  at  all  times  and  is  an 
input  variable  to  the  program.  The  non-isovelocity  environment  has  no  effect  on 
the  angular  horizontal  grid  cell  size  since  the  speed  of  sound  has  been  assumed 
to  be  a  function  of  depth  only.  In  general,  a  different-sized  patch  of  surface 
area  is  ensonified  by  each  surface  reverberation  cell.  Subroutine  TIMEBOUND 
is  usc'l  to  trace  the  upper  and  lower  bounding  rays  of  the  cell  to  the  surface, 
thus  determining  the  length  of  the  ensonified  patch.  The  width  of  the  patch, 
w,  is  determined  through  the  use  of  the  azimuthal  angular  increment,  Ad,  and 
the  radial  arm,  r,  from  the  vertical  axis  through  the  source  to  the  center  of  the 
scattering  element. 

w  =  rAd.  (3.8) 

Figure  3.19  displays  the  geometry  used  in  the  determination  of  the  grid  cell  width 
and  a  general  illustration  of  a  surface  reverberation  cell.  Once  the  ensonified 
area  (length  times  width)  of  the  surface  has  been  determined,  the  spectrum 
and  relative  power  gain  of  the  scattering  function  are  computed  through  the 
same  process  used  in  SCATTER.  The  subroutine  used  for  calculating  Doppler 
shift  in  SCATTER  has  been  modified  to  accept  different  angles  and  speeds  of 
sound  which  occur  at  different  scattering  elements.  The  new  subroutine  is  called 
RAYDOP  cind  a  subroutine  is  shown  in  Appendix  JI.13.  After  all  surfzu:e-solution 
sectors  have  been  processed,  the  volume  and  bottom  sectors  are  processed  in  a 
like  manner. 

Subroutine  VOLSEQ  is  used  to  determine  the  scattering  function  for  all 
volume  cells.  The  sectors  formed  from  the  volume  angles  are  each  divided  into 
cells  as  was  done  in  SURSEQ.  In  the  determination  of  the  backscattering  level 
for  each  volume  cell,  a  backscattering  coefficient  is  multiplied  by  the  volume  of 
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BOTTOM 


Figure  3.20  Vertical  wail  of  a  volume  reverberation  cell 
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each  cell.  This  volume  is  determined  by  tracing  rays  to  points  A,  B,  C,  and  D  as 
shown  in  Figure  3.20.  Subroutine  RAYSPOT  (Appendix  H.14),  which  utilizes 
subroutines  TURN  and  PIECES  (Appendices  H.15  and  H.16),  is  used  to  trace 
the  four  rays.  These  four  points  form  a  vertical  wall  of  the  cell,  which,  when 
multiplied  by  the  cell  width,  yields  the  volume  of  the  cell.  The  width  of  the 
volume  cells  is  determined  in  the  same  manner  as  the  width  of  the  surface  ceils, 
by  the  product  of  the  radial  arm,  r,  to  the  center  of  the  cell,  and  the  azimuthal 
incremental  angle,  A6.  An  algorithm  for  dividing  a  volume  cell  at  a  boundary, 
if  necessary,  is  also  implemented  in  VOLSEQ. 

In  the  final  SONAR  equation  subroutine,  BOTSEQ,  all  of  the  bottom- 
solution  sectors  are  processed  for  their  contributions  to  the  scattering  function. 
The  same  procedure  that  is  used  in  SURSEQ  is  followed  in  BOTSEQ,  with  one 
exception.  In  SURSEQ,  the  surface  wave  speed  is  used  to  obtain  an  additional 
Doppler  shift  component.  This  component  cannot  exist  for  the  bottom  of  the 
ocean  as  the  bottom  is  assumed  stationary.  Otherwise,  through  the  use  of 
the  bottom-solution  sectors  amd  ail  additional  bottom-related  parameters,  the 
contribution  by  the  ocean  bottom  to  the  scattering  function  is  determined  in  an 
algorithm  similar  to  that  in  SURSEQ. 

Once  the  scattering  function  for  each  type  of  reverberation  has  been 
calculated,  the  components  due  to  the  surface,  volume,  and  bottom  are  combined 
to  yield  the  total  scattering  function.  It  may  be  specified  beforehand,  however, 
that  only  one  type  of  reverberation  is  to  be  processed.  In  distinguishing  the 
different  types  of  reverberation,  the  convention  that  is  used  in  the  isovelocity 
model,  SCATTER,  is  maintained.  The  type  of  reverberation  is  dictated  by  the 
type  of  element  actually  doing  the  backscattering.  Surface  reverberation  has 
been  assumed  to  be  the  compilation  of  all  returns  which  are  backscattered  from 
the  surface,  even  though  a  certain  number  of  returns  may  have  reflected  from 
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the  bottom  several  times.  Likewise,  bottom  reverberation  is  comprised  of  all  of 
the  backscattered  energy  from  the  bottom.  All  of  the  remaining  reverberation 
is,  therefore,  due  to  volume  backscattering. 

3.2.3.  Capabilities  and  Output  from  RAYSCT 

The  revised  reverberation  modeling  program,  RAYSCT,  has  been  written 
to  perform  the  same  functions  as  SCATTER.  These  are  to  produce  either  a 
power  gain  versus  time  plot  or  a  spectral  density  plot  for  the  scattering  function. 
However,  many  more  physical  factors  have  been  taken  into  account  in  the  revised 
model  than  in  the  original.  A  sound  velocity  profile,  with  the  capability  of 
handling  up  to  fifty  constant  velocity  gradient  segments,  has  been  introduced. 
A  different  absorption  coefficient  is  represented  in  each  of  these  segments.  Also, 
multiple  forward  refiections  are  allowed,  along  with  associated  loss  factors  for 
surface  and  bottom  reflections.  These  additions  provide  a  simulation  of  more 
realistic  underwater  environments. 

The  output  from  RAYSCT  is  in  the  same  format  as  the  output  from 
SCATTER  as  described  in  section  3.1.3.  Therefore,  sample  plots  have  not  been 
included  in  this  section.  The  analysis  of  the  revised  model  has  been  reserved  for 
Chapter  4. 


Chapter  4 

COMPARISON  OF  REVERBERATION  MODELS 
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4.1.  Introduction 

The  results  presented  here  axe  comparisons  of  the  spectral  densities  and  power 
gains  of  scattering  functions  produced  by  the  original  model,  SCATTER,  and 
the  revised  model,  RAYSCT.  Because  of  the  very  large  number  of  parameters 
which  govern  the  output  from  these  models,  as  many  as  possible  have  been  kept 
constant,  limiting  the  number  of  comparisons  made.  There  are  many  features 
of  the  original  model  which  have  been  left  unchanged  in  the  revised  model; 
thus,  it  was  not  necessary  to  dwell  upon  them.  In  this  section,  the  presence  of 
multiple  boundary  returns  has  been  analyzed  with  respect  to  the  original  model. 
Additionally,  the  effects  of  different  sound  velocity  profiles  have  been  compared 
to  the  isovelocity  case. 

The  raypaths  for  all  surface  and  bottom  returns  in  this  chapter  are  referenced 
by  the  state  numbers,  KA  and  KB.  KA  has  been  defined  as  the  number  of 
times  a  ray  passes  through  the  region  between  the  source  and  the  surface,  and 
KB  Ba  the  number  of  times  a  ray  passes  through  the  region  between  the  source 
and  the  bottom.  The  state  numbers  that  correspond  to  a  particular  return 
define  the  raypath  to  that  boundary  from  the  source,  and  not  from  the  source 
to  the  boundary  and  back.  For  example,  if  a  bottom  return  has  occurred  after 
a  surface  reflection,  the  path  to  the  boundary  would  have  gone  from  the  source 
to  the  surface  and  then  to  the  bottom.  This  path  is  represented  by  the  state 
numbers,  KA  —  2,  and  KB  =  1,  and  is  denoted  by  “(2  —  1)”.  The  state  number 
corresponding  to  the  region  above  the  source,  KAy\a  always  listed  first. 
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4.2.  Description  of  the  Different  Scenarios  Used 

The  only  input  variables  to  the  new  model  that  differ  from  those  of  the 
original  model  are  the  surface  and  bottom  reflection  loss  factors.  For  the  scenarios 
chosen  as  examples,  four  parameters  were  varied.  These  include  the  source  depth, 
the  ocean  depth,  the  elapsed  time  between  transmitting  and  receiving,  and  the 
sound  velocity  profile.  Lists  of  the  input  variables  for  all  examples  to  be  shown 
appear  in  Appendix  E.  For  these  examples,  the  source  depth  may  be  either  one 
hundred  or  three  hundred  meters  below  the  surface,  and  the  ocean  depth,  either 
five  hundred  or  two  thousand  meters.  An  isovelocity  profile  is  used  for  some 
examples  and  a  sound  velocity  profile  consisting  of  both  a  negative  gradient 
layer  and  a  positive  gradient  layer  is  used  for  others. 

Two  types  of  environments  were  chosen  to  be  used  for  comparison.  The  first 
is  a  somewhat  shallow  ocean  with  an  isovelocity  profile.  It  will  serve  to  show 
the  results  of  the  inclxision  of  multiple  boundary  returns.  In  the  second  example, 
the  ocean  has  been  set  quite  a  bit  deeper  and  a  sound  velocity  profile  consisting 
of  a  decreasing  speed  of  sound  with  depth,  followed  by  an  increasing  speed  of 
sound  with  depth,  has  been  used.  Although  the  example  shows  some  effects  of 
multiple  boundary  returns,  the  main  emphasis  should  be  placed  on  the  effects  of 
the  sound  velocity  profile. 

4.3.  Discussion  of  Results 

4.3.1.  Multiple  Boundary  Returns 

A  scenario  which  includes  a  moving  source  has  been  shown  to  best  illustrate 
the  effects  of  the  surface  and  bottom  boundary  returns.  Therefore,  the  examples 
that  were  chosen  for  the  comparison  of  the  multiple  boundary  return  model  to 
the  single  boundary  return  model  include  a  moving  source. 
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The  first  scenario  used  for  comparison  is  identical  to  that  used  for  Figures  3.3 
and  3.5.  The  ocean  is  500  meters  deep  and  the  source  is  moving  with  a  velocity 
of  5  meters  per  second  at  a  depth  of  300  meters  below  the  surface.  The  speed  of 
sound  is  a  constant  1500  meters  per  second.  The  purpose  of  this  first  example  is 
to  compare  a  multiple  boundary  return  scattering  function  to  a  single  boundary 
return  scattering  function.  Raytracing  is  still  used,  but  refraction  is  not  a  factor. 

The  plots  generated  from  the  original  model,  SCATTER,  have  already  been 
discussed  in  Chapter  3.  Recall  the  power  gain  versus  time  plot  of  Figure  3.3  and 
the  single  spectrum  plot  of  Figure  3.5.  In  the  single  spectrum  plot,  the  surface 
and  bottom  contributions  aopear  at  different  Doppler  shift  values  implying  a 
different  angle  of  backscattering  from  the  surface  than  from  the  bottom.  The 
degree  of  Doppler  shifting  of  the  surface  and  bottom  returns  will  also  be  shown 
to  be  very  important  in  the  analysis  of  the  multiple  boundary  return  examples. 

For  the  comparable  RAYSCT  example,  all  input  parameters  were  kept  the 
same  as  in  the  original  model.  The  only  input  variables  added  to  the  revised 
model  are  the  loss  factors  for  surface  and  bottom  reflections.  For  this  example, 
ray  tracing  is  used  to  demonstrate  multiple  boundary  returns,  but  as  the  velocity 
gradient  has  been  set  equal  to  zero,  straight-line  propagation  results.  The  input 
parameters  for  the  first  RAYSCT  scenario  may  be  found  in  Appendix  E.Z. 

Figure  4.1  contains  the  power  gain  versus  time  ciirves  for  the  surface,  volume, 
bottom,  and  combined  total  contributions  to  the  scattering  function.  The 
additional  surface  and  bottom  boundary  returns  are  seen  to  cause  a  slightly 
slower  decay  in  the  power  gain  curve  as  compared  to  the  power  gain  curve  of 
Figure  3.3.  The  individual  surface  and  bottom  contributions  also  decay  more 
slowly  than  in  Figure  3.3  because  of  the  increased  energy  from  the  additional 
boundary  returns.  The  volume  reverberation,  also,  has  a  lower  decay  rate 
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associated  with  it.  This  is  due  to  the  additional  volume  reverberation  originating 
after  the  occurrence  of  boundary  returns. 

It  might  be  observed  that  the  first  occurrences  of  surface  and  bottom  returns 
in  Figure  4.1  appear  at  a  slightly  lower  level  than  those  in  Figure  3.3.  This 
is  because  the  approximations  for  the  ensonified  areas  of  surface  and  bottom 
cells  is  more  accurate  in  the  revised  model,  RAYSCT.  The  width  of  the  cell  is 
approximated  by  the  arc  across  the  middle  of  the  grid  ceil  in  RAYSCT  (Figure 
3.22),  as  opposed  to  the  arc  at  the  front  (smaller  end)  of  the  grid  cell  as  is  assumed 
in  SCATTER.  The  approximation  in  RAYSCT  yields  larger,  but  more  accurately 
represented  grid  cells,  which  directly  accounts  for  a  slightly  lower  intensity  than 
the  approximation  in  SCATTER. 

Figure  4.2  is  a  spectral  density  plot  taken  at  a  time  of  1.74  seconds.  The 
peak  at  zero  Doppler  shift  is  due  to  a  combination  of  the  first  surface  (1  —  O)  and 
the  first  bottom  (0  —  1)  returns.  These  returns  occur  at  or  near  zero  Doppler 
because  the  time  at  which  the  spectrum  was  taken  is  such  that  the  grazing  angles 
of  these  returns  at  the  boundaries  are  quite  low.  The  next  peak  to  the  left,  at 
approximately  —60  Hz  Doppler  shift,  is  due  to  the  next  pair  of  surface  (1  —  2) 
and  bottom  (2  —  l)  returns.  These  returns  are  just  approaching  the  zero  Doppler 
point  and  appear  in  the  same  frequency  region  as  each  other.  The  next  frequency 
peak  in  Figure  4.2,  occurring  at  about  —170  Hz  Doppler  shift,  is  solely  the  third 
bottom  (2  —  3)  return.  The  third  surface  (3  —  2)  return  with  which  it  is  paired, 
has  not  quite  merged  with  it  to  the  same  frequency  region.  At  this  time  (t  =  1.74 
seconds),  the  volume  contribution  may  be  observed  at  low  levels  between  —100 
Hz  and  +30  Hz  Doppler  shift. 
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4.3,2,  Non-isovelocity  Profile 

Raytracing  may  have  drastic  effects  on  the  scattering  function.  The  scenario 
in  the  next  example  presented  consists  of  a  longer  time  interval  than  in  the  first 
set  of  examples.  The  source  is  situated  one  hundred  meters  from  the  surface  at  a 
2000-meter  ocean  depth.  The  source  has  a  velocity  of  five  meters  per  second  and 
ail  scattering  elements  have  a  mean  velocity  of  zero.  One  set  of  plots  has  been 
generated  by  SCATTER.  In  this  set,  an  isovelocity  profile  and  no  reflections 
at  the  boundaries  have  been  assumed.  A  comparative  set  of  plots  has  been 
generated  by  RAYSCT,  the  modified  program.  For  these  plots,  the  options  of  a 
non-isovelocity  profile  and  the  allowance  of  multiple  boundary  returns  have  been 
exercised. 

A  power  gain  versus  time  plot  is  shown  in  Figure  4.3  for  the  surface,  volume, 
bottom,  and  total  contributions  to  the  scattering  function  for  the  isovelocity  case, 
SCATTER.  The  input  variables  that  were  used  are  shown  in  Appendix  EA.  The 
curves  in  Figure  4.3  are  quite  similar  to  the  curves  which  appeared  in  Figure  3.3. 
The  major  difference  is  that,  in  Figure  4.3,  the  surface  return,  which  has  occurred 
almost  immediately  after  transmitting,  now  precedes  the  bottom  return.  This 
bottom  return,  it  can  be  seen  from  the  figure,  appears  at  a  time  equal  to  about 
2.5  seconds.  The  times  corresponding  to  the  surface  and  bottom  returns  result 
directly  from  the  situation  of  the  source  in  the  ocean;  it  is  very  close  to  the 
surface,  and  relatively  far  from  the  bottom.  The  surface  reverberation  seems  to 
dominate  for  the  entire  range  shown. 

For  the  non-isovelocity  case,  the  scattering  function  proves  to  be  quite 
different.  The  sound  velocity  profile  that  was  used  is  shown  in  Figure  4.4.  This 
type  of  profile  might  be  the  result  of  a  warm  ocean  surface  suspended  above  a 
colder  layer  of  water.  Thus,  the  decrease  in  temperature,  and  also  sound  speed, 
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is  evident  to  a  depth  of  two  hundred  meters.  The  increase  in  the  speed  of  sound 
with  depth  below  two  hundred  meters  would  be  due  to  the  pressure  increase  with 
depth. 

Power  gain  versus  time  plots  for  the  surface,  volume,  bottom,  and  total 
contributions  to  the  scattering  function  are  shown  in  Figure  4.5.  The  input 
parameters  that  were  used  are  shown  in  Appendix  E.5.  The  surface  reverberation 
has  actually  disappeared  after  a  time  of  about  1.3  seconds.  This  occurs  during 
the  (1  —  0)  surface  return.  The  ray  diagram  shown  in  Figure  4.6  explains  this 
phenomenon.  The  rays  have  been  traced  for  2.5  seconds  in  two  clusters  of 
equispaced  transmit  angles:  9°  <  ^  <  13°  in  one  half  degree  increments  and 
—88°  <  4>  <  —80°  in  two  degree  increments.  At  the  time  at  which  the  surface 
return  vanishes,  there  is  no  possible  ray  solution  to  the  surface;  thus,  there  can 
be  no  surface  reverberation.  A  critical  ray,  transmitted  at  11.3  degrees,  grazes 
the  surface  at  a  time  of  0.66  seconds.  Any  ray  transmitted  at  an  angle  higher  in 
elevation  than  this  critical  angle,  is  reflected  from  the  surface  and  yields  volume 
reverberation.  Any  ray  transmitted  below  this  critical  angle  is  refracted  away 
from  the  surface  and  also  yields  volume  reverberation.  This  refraction  of  energy 
away  from  the  surface  forms  a  shadow  zone,  a  region  in  which  the  acoustic 
intensity  is  extremely  weak  compared  to  neighboring  regions.  The  shadow  zone, 
which  is  shown  in  Figure  4.6,  begins  at  a  time  of  0.66  seconds.  This  shadow  zone 
directly  accounts  for  the  loss  in  surface  reverberation  at  a  time  of  1.3  seconds. 
Although  the  acoustic  intensity  is  weak  in  a  shadow  zone.  Figure  4.5  shows  the 
model  predicts  that  the  surface  retxirn  will  vanish  completely.  This  is  because 
ray  theory  predicts  zero  intensity  in  a  shadow  zone. 

Figure  4.5  also  shows  that  the  surface  reverberation  reappears  at  some  later 
time.  The  surface  return  is  rejuvenated,  not  through  a  direct  path  to  the  surface, 
but  from  a  bottom  reflection.  Although  the  direct  surface  return  has  vanished. 


The  ray  diagram  illustrates  the  emergence  of  a  second 
surface  return  through  a  bottom  reflection. 
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the  somewhat  attenuated  (1  —  2)  source  to  bottom  to  surface  return  is  very 
evident  at  a  time  of  about  5.1  seconds.  This  surface  return  is  weak  since  it  has 
undergone  two  bottom  reflections,  one  on  the  way  to  the  scatterer,  and  one  along 
the  same  path  back  to  the  source.  In  fact,  the  reverberation  is  dominated  by  the 
bottom  contribution  at  this  point. 

Another  peculiarity  of  Figure  4.5  may  be  observed  in  the  volume  component 
of  the  scattering  function  beyond  the  approximate  4.9  second  mark.  The  power 
gain  has  increased  rapidly  and  the  curve  has  lost  its  smooth,  continuous  decay. 
For  investigative  purposes,  a  ray  diagram  of  a  portion  of  the  volume  has  been 
generated  as  shown  in  Figure  4.7.  The  rays  in  the  range  of  —10®  <  (f>  <  +10® 
have  been  traced  for  8  seconds.  It  can  be  seen  that  a  strong  caustic  has  formed 
at  a  distance  of  about  4000  meters.  This  would  correspond  to  an  approximate 
travel  time  of  5.3  seconds  to  the  volume  scattering  elements  and  back  to  the 
receiver.  The  catistic  region  actually  began  at  a  distance  of  about  3500  meters 
and  extends  well  beyond  this.  The  volume  contribution  to  the  scattering  function 
reflects  this  caustic.  Additionally,  the  rays  involved  in  the  formation  of  the  caustic 
have  refracted  away  from  the  surface  and  the  bottom,  and  therefore,  have  not 
experienced  any  reflection  loss  to  these  boundaries. 

The  increase  in  power  gain  due  to  the  caustic  may  be  explained  as  follows. 
The  rays  which  vertically  bovmd  the  cells  involved  in  the  caustic  are  shown  to  be 
converging  in  Figure  4.7,  accounting  for  a  very  high  intensity  level  at  these  cells. 
But,  because  the  rays  are  very  closely  spaced,  the  volumes  of  these  cells  are  small, 
accounting  for  relatively  low  backscattered  levels.  These  low  backscattered  levels 
act  to  negate  the  high  intensity  levels  incident  on  the  volume  scattering  elements. 
The  increased  power  gain  shown  in  Figure  4.5  has,  to  this  point,  not  yet  been 
explained.  It  may  not  simply  be  attributed  to  a  one-way  focusing  of  energy  on 
some  localized  region  out  in  the  field.  The  principle  of  acoustic  reciprocity  may 
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be  used  to  shed  some  light  on  the  problem.  It  states  that,  “if  in  an  unchanging 
environment,  the  locations  of  a  small  source  and  a  small  receiver  are  interchanged, 
the  received  signal  will  remain  the  same”  [6].  Utilizing  this  principle,  it  becomes 
clear  that  the  backscattered  sound  from  each  cell  involved  in  the  caustic  is  focused 
back  on  the  receiver  in  the  same  degree  that  the  transmitted  sound  had  been 
focused  on  the  scatterer.  It  is  this  additional  focusing  of  sound  on  the  return 
path  by  many  cells  which  is  the  cause  for  higher  power  gain  levels  at  the  receiver. 
Therefore,  the  rise  in  power  gain  that  is  shown  in  Figure  4.5  as  the  caustic  is 
approached  in  time,  should  be  expected. 
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Chapter  5 

SUMMARY,  CONCLUSIONS,  AND  RECOMMENDATIONS 
5.1.  Summary 

Multiple  boundary  retxirns  and  non-isovelocity  gradients  can  have  significant 
effects  on  underwater  acoustic  reverberation  originating  from  a  moving  pulsed 
source.  An  isovelocity  reverberation  model,  which  disallows  surface  and  bottom 
boundary  reflections  has  been  modified  to  include  raytracing  and  multiple  bound¬ 
ary  returns.  The  portion  of  the  model  which  was  revised  is  that  used  to  compute 
the  ocean  scattering  function  which  characterizes  the  received  reverberation.  It 
is  independent  of  the  transmitted  pulse  shape  as  well  as  its  power  level.  Effects 
from  the  source  and  scatterer  motion,  the  backscattering  cross-section  of  the 
ocean  elements,  the  acoustic  beam  effects,  aoid  the  propagation  losses  are  all 
considered  in  the  determination  of  the  scattering  function. 

In  the  modified  scattering  function  computation,  allowances  for  reflections 
from  the  air-water  interface  and  the  ocean  bottom  have  been  made.  Reflection 
loss  factors  for  surface  and  bottom  reflections  have  been  made  input  parameters 
to  the  revised  model.  Additionally,  a  depth-dependent  sound  velocity  profile  has 
been  incorporated.  As  many  as  fifty  depth-velocity  pairs  of  points  are  allowed. 
The  sound  velocity  profile  has  been  assumed  to  consist  of  constant  velocity 
gradient  segments.  The  consequence  of  using  such  straight  line  segments  is  that 
the  raypaths  through  any  one  horizontal  stratification  are  arcs  of  circles. 

A  series  of  raytracing  routines  has  been  written  to  perform  all  raytracing 
functions  required  for  the  modification  of  the  original  scattering  function 
program.  The  raytracing  theory  that  has  been  used  for  the  present  model  has 
been  based  on  the  theory  from  Clay  and  Medwin  [2]. 
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The  raytracing  routines  have  been  written  to  perform  two  major  functions. 
First,  given  a  travel  time,  t,  and  an  initial  transmitted  angle,  0^,  the  final  position 
of  the  ray  is  determined.  The  significance  of  this  is  that  a  constant  time  contour, 
over  which  to  build  a  grid  of  cells,  may  be  defined.  In  the  original  isovelocity 
model,  the  scattering  functions  from  a  spherical  grid  of  cells  are  summed  to 
obtain  the  total  scattering  function.  A  sphere  is  actually  a  constant  distance 
contour.  However,  since  the  original  model  has  assumed  a  constant  sound  speed, 
the  sphere  is  also  a  constant  time  contour.  When  a  non-isovelocity  environment 
is  introduced,  the  constant  time  contour  no  longer  takes  the  shape  of  a  sphere. 
The  actual  constant  time  contour  is  defined  by  the  first  raytracing  function 
which  yields  the  final  position  of  a  ray  given  the  initial  angle  and  the  travel 
time.  This  raytracing  routine  also  determines  the  final  angle,  path  length, 
and  absorption  loss  of  the  ray,  as  well  as  the  numbers  of  surface  and  bottom 
reflections.  Information  from  this  routine  is  used  to  determine  the  reflection  loss 
and  the  backscattering  cross-section  of  any  grid  cell. 

The  second  purpose  of  the  raytracing  package  is  to  define  all  surface  and 
bottom  reverberation  cells  at  a  particular  time.  The  original  model  has  zussumed 
that  the  only  surface  reverberation  cells  are  those  on  the  direct  path  from  the 
source  to  the  surface,  and  that  the  only  bottom  returns  are  on  the  direct  path  to 
the  bottom.  However,  if  the  travel  time  is  long  enough,  there  may  be  a  number 
of  groups  of  both  surface  and  bottom  reverberation  cells.  In  the  raytracing 
routine  in  which  the  surface  and  bottom  cells  are  defined,  the  travel  time  and  the 
boundary  of  interest,  either  the  surface  or  the  bottom,  are  the  input  parameters. 
The  angles,  of  all  rays  which  end  at  the  boundary  at  the  specified  time  are 
output  from  the  routine.  When  the  pulse  length  is  known,  this  routine  may  be 
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used  to  determine  all  groups  of  transmit  angles  at  the  source  which  eventually 
yield  surface  or  bottom  reverberation  cells.  This  routine  is  also  used  to  determine 
the  areas  of  the  ensonified  surface  and  bottom  regions  at  a  particular  time. 

The  scattering  function  output  from  the  original  model  assumed  the  form 
of  either  the  power  gain  versus  time  or  the  spectral  density  of  the  scattering 
function  for  one  particular  time.  The  same  convention  is  used  for  the  output 
from  the  revised  model,  RAYSCT. 

5.2.  Conclusions 

The  inclusion  of  a  non-isovelocity  profile  and  multiple  boundary  returns  to 
an  isovelocity  reverberation  model  has  been  shown  to  add  a  significant  amoimt 
of  detail  to  the  scattering  function  of  the  ocean.  Because  of  the  allowance  of 
reflections,  the  decay  of  the  volume  reverberation  was  not  as  rapid  in  the  modified 
model.  The  additional  surface  and  bottom  returns  have  been  shown  in  many 
cases  to  maintain  the  overall  reverberation  at  a  higher  level  in  the  revised  model. 

Even  more  pronounced  differences  have  been  shown  in  the  frequency  spectra 
of  the  scattering  functions.  When  a  sotirce  is  in  motion,  a  frequency  offset 
is  introduced  for  all  scatterers  that  are  off  the  axis  of  source  motion.  Thus, 
whenever  a  boundary  return  is  observed  for  the  first  time  on  a  scattering  function 
versus  time  plot,  it  is  maximally  Doppler  shifted,  as  its  direction  of  arrival  is 
ninety  degrees  off-axis  of  the  source  motion.  For  the  cases  that  were  presented, 
the  multiple  boundary  returns  were  not  shown  to  contribute  enormously  to  the 
power  gain  curves,  although,  some  significant  differences  were  noted.  However, 
the  multiple  boundary  returns  have  been  shown  to  contribute  very  significantly 
to  the  scattering  function  spectra. 

When  a  non-isovelocity  environment  was  introduced,  the  power  gain  curves, 
as  well  as  the  spectral  density  plots  were  sometimes  observed  to  be  extremely 
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different  from  the  plots  generated  by  the  isovelocity  model.  In  the  example 
chosen,  the  surface  reverberation  went  away  completely  after  a  certain  length  of 
time.  This  occurrence  was  shown,  via  raytracing  diagrams,  to  be  the  result  of  a 
shadow  zone  at  the  surface.  The  surface  component  of  the  scattering  function 
reappeared  at  a  later  time  due  to  a  different  surface  return.  The  path  to  the 
surface  for  the  second  return  included  a  bottom  reflection. 

The  particular  example  that  was  iised  to  demonstrate  ray  effects  did  not  show 
drastically  different  results  than  were  expected  for  the  bottom  contribution  to  the 
scattering  function.  However,  because  the  surface  return  had  disappeared  and 
then  had  returned  at  a  lower  level,  the  bottom  component  was  seen  to  dominate 
the  power  gain  versus  time  curve,  even  after  the  surface  component  had  returned. 
This  is  an  expected  result  because,  although  the  surface  backscattering  coefficient 
had  been  input  as  a  larger  value  than  the  bottom  backscattering  coefficient,  the 
surface  component  had  been  attenuated  during  bottom  reflection. 

The  volume  reverberation  from  the  non-isovelocity  scenario  has  shown  some 
very  pronounced  effects  from  a  caustic  region.  The  volume  scattering  function 
spectrum  versus  time  was  seen  to  increase  at  a  time  corresponding  to  the 
occurrence  of  the  caustic.  A  ray  diagram  has  been  used  to  conflrm  the  existence 
of  this  caustic.  The  increase  in  power  gain  was  attributed  to  the  focusing 
of  backscattered  energy  from  volume  elements  onto  the  receiver  at  the  source 
location.  The  principle  of  acoustic  reciprocity  was  used  to  point  out  that  the 
focusing  of  energy  is  the  same  from  a  scattering  element  baM:k  to  the  receiver  as 
it  had  been  from  the  source  to  that  scattering  element.  Therefore,  the  power 
gain  in  the  case  of  a  caustic  was  shown  to  be  justifiably  different  from  the  power 
gain  in  an  isovelocity  scenario,  in  which  spherical  spreading  exists  along  all  paths 
both  to  and  from  all  elements. 
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5.3.  Recommendations  for  Further  Work 

In  this  effort,  ray  tracing  and  multiple  boundary  returns  have  been  introduced 
to  an  isovelocity  reverberation  model.  Specular  reflections,  those  in  which  the 
angle  of  reflection  is  equal  to  the  angle  of  incidence,  have  been  assumed.  All 
reverberation  has  been  assumed  to  result  from  backscattered  energy  only;  the 
path  from  the  scatterer  back  to  the  receiver  is  the  same  path  as  that  from  the 
source  to  the  scatterer.  This  is  referred  to  as  the  monostatic  case.  There  are 
additional  source  to  receiver  paths  which  may  exist,  but  they  had  been  considered 
to  be  out  of  the  scope  of  this  effort.  These  reverberation  paths,  in  which  the 
return  path  is  not  the  same  as  the  path  to  the  scattering  element,  may  be  referred 
to  as  bistatic  returns.  An  example  of  such  a  return  would  be  the  raypath  from 
the  source  to  the  surface,  then  to  the  bottom  after  a  surface  reflection,  and  finally 
back  to  the  source  from  the  bottom.  This  class  of  returns  should  be  considered 
in  the  future.  Some  of  the  parameters  of  the  model  that  would  have  to  be 
examined  if  bistatic  returns  were  to  be  incorporated  are  the  Doppler  shift,  the 
grid  cell  depth,  the  bistatic  scatterer  radiation  pattern  and  the  acoustic  beam 
effects.  The  Doppler  due  to  a  moving  source  would  involve  two  factors:  one 
resulting  from  the  transmit  angle  and  one  resulting  from  the  receive  angle.  The 
grid  cell  depth  could  no  longer  be  defined  as  simply  r/2,  in  seconds,  because  the 
energy  would  not  be,  in  general,  directly  backscattered,  but  would  be  scattered  at 
some  other  angle.  The  determination  of  the  backscattering  cross  section  should 
follow  directly  from  the  grid  cell  depth.  Finally,  the  beam  weighting  factor  would 
have  to  consist  of  a  transmit  beam  attenuation  at  one  angle,  and  a  receive  beam 
attenuation  at  some  other  angle. 

Other  work  which  would  be  beneficial  to  the  model,  would  be  the  incor¬ 
poration  of  frequency  and  grazing  angle  dependent  scattering  coefficients  and 
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reflection  loss  factors.  The  FORTRAN  computer  code,  which  was  written  in  this 
present  effort,  has  been  left  open  so  that  such  modifications  can  be  made. 

Finally,  the  model  should  be  tested  for  validity  by  running  some  comparisons 
to  in-water  data.  For  this  to  be  done,  actual  acoustic  beaun  patterns  would  have 
to  be  \ised,  a  transmit  signal  envelope  and  power  level  would  have  to  be  known, 
as  well  as  some  detailed  information  on  the  particular  environment  used.  This 
would  include  an  accurate  sound  velocity  profile,  the  wind  speed,  the  bottom 
type,  the  sea  state,  and  the  existence  of  any  moving  current  layers.  The  revised 
model’s  performance  on  the  prediction  of  reverberation  from  caustics  and  shadow 
zones  should  also  be  verified. 
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APPENDIX  A 

DOPPLER  SHIFT  CALCULATIONS 
A.l.  Source  and  Receiver  Motion 

Let  V  represent  the  source  velocity.  The  Doppler,  resulting  at  an  angle,  7, 
from  the  axis  of  source  motion,  is  related  to  the  component  of  the  source  velocity 
in  that  direction.  This  component,  v,  is  equal  to  V  cos  7. 

During  a  pulse  of  duration,  r  seconds,  and  frequency,  /o,  there  exist  N  cycles. 
The  following  relation  holds: 

.  (A.1) 

T 

The  length  of  the  pulse  is  cr  in  meters,  where  c  is  the  speed  of  sound  in  meters 
per  second.  Consider  the  leading  edge  of  a  pulse,  transmitted  at  a  time,  to,  from 
a  moving  sottrce,  at  some  point,  Zo.  By  the  time  the  trailing  edge  is  transmitted, 
the  leading  edge  will  have  traveled  a  distance  of  ct  from  its  point  of  origin,  Z0. 
The  source,  itself,  however,  will  have  traveled  a  distance  of  vr,  where  v  =  V  cos  7. 
The  length  of  the  transmitted  pulse  in  the  water  is  the  difference  between  the 
leading  and  trailing  edges  of  the  pulse. 


CT  =  cr  —  wr  =  (c  —  v)t 


(A.2) 


where  is  the  duration  of  the  Doppler-shifted,  transmitted  pulse, 
equation  A.2, 


_  (c  -  v)r 

D-  c  - 

The  Doppler-shifted  frequency  may  be  represented  as 


From 

(A.3) 


/ 

J  n  -  • 


Substitution  of  equation  A.3  into  A.4  yields 

(r^) 


(A.4) 


(A.5) 
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This  may  be  rewritten  as 

=  fo  (r^)  ^  c  ^  ^ 

For  V  c,  the  higher  ordered  terms  of  the  series  in  equation  A.6  may  be 
neglected,  yielding 

/^  =  fo  +  ^cos'y^  .  (A.7) 

For  the  case  of  a  moving  receiver,  the  same  reasoning  applies.  If  returns 
to  the  receiver  at  the  same  angle,  7,  off  of  the  axis  of  receiver  motion,  then  the 
received  frequency,  /^,  may  be  expressed  as 

Again,  if  the  higher  ordered  terms  are  considered  negligible,  the  received 
frequency  becomes 

^cos7^  .  (A.9) 

A.2.  Scatterer  Motion 

A  somewhat  different  situation  arises  when  the  transmitted  signal  is  backscat- 
tered  from  some  moving  point,  z®.  The  pulse,  of  frequency,  fp,  begins  with  a 
length  of  cr  before  backscattering.  Consider  a  scatterer,  moving  with  a  velocity, 
U,  at  some  angle,  7,  from  the  direction  of  the  approzudiing  pulse.  The  component 
of  the  scatterer  velocity  in  the  direction  of  the  pulse,  u,  is  equal  to  17  cos  7. 

At  some  time,  to,  the  leading  edge  of  the  pube  reaches  the  scattering  point, 
Xg.  As  the  pulse  folds  back  over  itself,  the  trailing  edge  reaches  the  new  scatterer 
position,  j/g,  in  some  time  At,  after  the  front  of  the  pulse  had  reached  the 
scatterer.  The  distance  traveled  by  the  trailing  edge  in  the  time.  At,  is 


Azi  =  cAt. 


(A.IO) 
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See  Figure  A.l.  The  leading  edge  of  the  pulse  is  labeled,  L,  and  the  trailing  edge, 
T.  In  the  same  amount  of  time.  At  the  scatterer  has  moved  from  Xg  to  a 
distance  of 

Axj  =  uAt.  (-^*11) 

From  the  figure,  it  can  be  seen  that  the  original  pulse  length,  cr,  is  equal  to  the 
sum  of  Axi  and  Ax2. 

cr  =  cAt  +  uAt  =  (c  +  u)  At.  (■^•12) 

The  length  after  backscattering,  cr^,  is  the  difference  between  the  leading 
and  trailing  edges  of  the  pulse  after  At  has  passed.  The  distance  of  the  leading 
edge  of  the  pulse  from  the  original  point,  x^,  at  time.  At,  is 

XI  =  cAt,  (A.13) 

and  the  distance  of  the  trailing  edge  from  the  same  point,  is 

xj-  =  uAt.  (■^•14) 

The  difference  in  these  distances  is  the  Doppler-shifted  pulse  length. 

CT  =  cAt  —  uAt  =  (c  —  u)At.  (A.15) 

From  equation  A.  12, 

At  =  (A.16) 

(c  +  tt) 

Whereas,  from  equation  A.  15, 

At  =  {A.17) 

(c-u) 

Equating  equations  A.16  and  A.17, 

,  _  ~  “) 

D-  (c  +  u)  • 


(A.18) 
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Incorporating  equations  A.l  and  AA  into  A.18,  the  backscattered,  Doppler- 
shifted  frequency  may  be  expressed  as 


(A.19) 


This  may  be  rewritten  as 


Finally,  if  u  <C  c,  the  higher  ordered  terms  may  be  neglected  and  the  following 
result  obtained: 

^  ^  •  (^-21) 


I 
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APPENDIX  C 

LIST  OF  INPUT  VARIABLES 

V  S  source  speed 

V  D  source  depth 

DOC E AN  ocean  depth 

ALPHA  absorption  coefficient 
SS  sound  speed 

SSCAT  surface  scattering  coefficient 
SSIG  standard  deviation  of  surface  scatterer  speed  distribution 
SWT  HE  surface  wave  vector  of  motion 

SW S  surface  wave  speed 

SWSIG  standard  deviation  of  stirface  wave  speed  distribution 

RATIO  power  ratio  of  surface  wave  to  surface  scattering  component 
VSCAT  volume  scattering  coefficient 

VSIG  standard  deviation  of  volume  scatterer  speed  distribution 
BSC  AT  bottom  scattering  coefficient 

BSIG  standard  deviation  of  bottom  scatterer  speed  distribution 
NCL  number  of  current  layers 

CTHICK{I)  thickness  of  current  layer,  I 
CTHE{I)  direction  of  motion  of  current  layer,  I 
CS{I)  speed  of  current  layer,  I 
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FREQ  carrier  frequency 
TPLL  pulse  length 

TTHE  transmit  beam  vector  in  azimuth 

TPHI  transmit  beam  vector  in  elevation 

RTHE  receive  beam  vector  in  azimuth 
RPHI  receive  beam  vector  in  elevation 

BRING  beam  pattern  file  angular  increment 
NBP  number  of  beam  pattern  points 

FS  baseband  sampling  rate 

NFFT  number  of  bins  across  bandwidth 

BELT  HE  grid  cell  width  in  azimuth 
TSTHE  total  spread  of  grid  azimuth  centered  on  receive  beam  axis 
DELPHI  grid  cell  height  in  elevation 
TSPHI  total  spread  of  grid  elevation  centered  on  receive  beam  axis 
RB  range  beginning 

RE  ramge  ending 

RI  range  increment 

OUTCNl  scattered  spectrum  or  total  power  flag 


OUTCN2  surface,  volume,  bottom,  total,  or  all  four  flag 
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APPENDIX  D 

ISOVELOCITY  MODEL  SUBROUTINES 
D.l.  Initialization  Subroutines 
D.1.1.  INRVMD 

Subroutine  INRVMD  simply  reads  into  SCATTER  the  input  parameters 
from  the  direct  access  file  created  back  in  INPUT.  If  it  is  desired,  these  input 
parameters  may  be  displayed.  There  is  also  the  option  in  this  subroutine  to  use 
two-dimensional  beam  pattern  weightings,  as  opposed  to  the  assumed  circularly 
symmetric  ones  created  in  BEAMS. 

D.1.2.  INITRM 

Subroutine  INITRM  prepares  the  input  parameters  for  use  in  program 
SCATTER.  It  is  called  once  each  time  SCATTER  is  run.  In  this  subroutine, 
the  backscattering  coefiicients  are  converted  from  decibels  to  linear  units;  the 
standard  deviation  of  scatterer  speed  distributions  are  converted  from  meters 
per  second  to  hertz;  and  all  angles  are  converted  from  degrees  to  radians.  The 
variables  defining  the  grid  are  initialized  zind  the  arrays  for  the  scattering  function 
accumulations  are  initialized  to  zero.  Also,  the  output  file  for  use  in  TRANSMIT 
is  defined,  and  a  header,  a  list  of  constants,  is  written  to  this  output  file. 

D.2.  Determination  of  Surface,  Volume,  or  Bottom  Reverberation 

D.2.1.  VANGLE 

Subroutine  VANGLE  is  called  once  for  each  range.  Here,  the  vertical  angles 
from  the  source  to  the  surface  and  to  the  bottom  at  both  the  front  and  the 
back  of  the  grid  cell  depth  are  calculated.  The  front  of  the  cell  is  defined  as  the 
distance  corresponding  to  the  range  and  the  back  of  the  cell  corresponds  to  the 
range  plus  the  grid  cell  depth.  Since  this  model  does  not  consider  any  refiections 
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from  the  surface  or  the  bottom,  the  first  vertical  cell  in  the  grid  is  at  the  surface 
where  the  range  (in  meters)  reaches  the  surface.  If  the  current  range  is  such  that 
the  surface  is  not  yet  reached  by  the  pulse,  then  the  first  vertical  cell  appears 
straight  upward  from  the  source,  or  at  <^  =  +90°.  A  similar  situation  exists  for 
the  bottom  with  the  last  vertical  cell.  If  the  bottom  is  farther  from  the  source 
than  the  current  range,  then  the  last  vertical  cell  appears  straight  downward  at 
<i>  =  —90°.  Subroutine  VANGLE  utilizes  subroutine  TOTANG  to  determine  the 
total  number  of  vertical  angular  increments,  or  rows,  in  the  grid.  First  found 
aire  the  numbers  of  rows  of  cells  pertaining  to  each,  surface,  volume,  and  bottom 
reverberation.  The  sum  of  these  three  is  equal  to  the  total  number  of  rows  in 
the  grid. 

D.2.2.  INITDO 

Program  SCATTER  has  the  option  of  computing  either  surface,  volume, 
bottom,  or  total  reverberation.  It  may  also  compute  all  four.  It  will  have  been 
predetermined  (in  program  INPUT)  exactly  which  type  of  reverberation  is  being 
sought.  Subroutine  INITDO  sets  the  “do  loop”  parameters  for  the  number  of 
rows  of  cells  to  be  considered.  For  example,  if  only  surface  reverberation  is 
of  Interest,  then  the  reverberant  power  gain  from  only  those  rows  of  grid  cells 
which  fall  on  the  surface  is  to  be  accumulated.  If,  on  the  other  hand,  the  total 
reverberation  spectrum  is  desired,  then  the  power  gain  from  all  rows  of  cells  is  to 
be  accumulated.  Subroutine  INITDO  prescribes  the  first  and  last  rows  of  cells 
over  which  to  accumulate  the  reverberation  power  gain. 

D.3.  Beam  Patterns  -  Subroutine  BPATRM 

The  total  reverberation  level  from  each  cell  is  weighted  by  both  the 
transmitting  and  receiving  beams.  As  was  stated  earlier,  these  beams  may  each 
be  steered  in  any  direction  through  the  specification  of  a  steering  vector,  {0,4)y 
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where  0  is  the  azimuthal  angle  and  ^  is  the  elevation  angle.  Subroutine  BPATRM 
computes  the  beam  weightings  for  each  ceil  by  interpolating  between  the  beam 
pattern  points  that  had  been  input.  Adjustments  are  made  for  the  steering 
vectors  for  each  beam  and  the  transmit  and  receive  beam  pattern  weightings  are 
multiplied  together  to  obtain  a  composite  weighting  at  the  receiver  for  that  cell. 

D.4.  Frequency  Spreading  and  Shifting  Routines 

D.4.1.  ZPDF 

Program  SCATTER  accumulates  the  reverberation  power  gain  by  beginning 
with  a  row  of  cells,  computing  the  sonar  equation  for  each  cell  in  that  row,  and 
then  proceeding  to  the  next  row  to  be  processed.  As  it  is  structured,  SCATTER 
steps  through  all  surface  reverberation  cells  first,  then  volume,  then  bottom. 
Subroutine  ZPDF  is  called  only  for  the  first  row  of  cells  involved  in  each  type  of 
reverberation.  Thus,  it  is  called  once  for  the  surface,  once  for  the  volume,  and 
once  for  the  bottom. 

The  actual  function  of  subroutine  ZPDF  is  to  calculate  the  quantized  and 
truncated  zero  mean  Gaiissian  scatterer  velocity  distribution.  The  probability 
density  function  defining  the  random  scatterer  motion  is  assumed  to  be  a 
Gaussian  curve  shaped  by  the  standsu'd  deviation,  a.  The  standard  deviations  of 
the  surface,  volume,  and  bottom  scatterer  speed  distributions,  which  had  been 
converted  from  m/s  to  Hz  are  input  variables  to  INPUT.  These  parameters 
are  what  determine  the  amount  of  frequency  spreading  which  occurs  at  each 
backscattering  element.  A  standard  deviation  of  0.0  would  denote  that  the 
scattering  element  has  no  motion  at  all,  and  that  the  backscattered  spectral 
distribution  would  be  identical  to  the  incident  spectral  distribution. 
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D.4.2.  DOPP 

Subroutine  DOPP  calculates  two  Doppler  shift  components  for  each  ceil  and 
adds  a  third.  First  of  all,  the  Doppler  due  to  source  motion  relative  to  the 
cell  is  determined.  Secondly,  the  Doppler  attributed  to  the  mean  scatterer  speed 
relative  to  a  stationary  source  is  found.  Additionally,  a  frequency  offset,  such  that 
a  stationary  scatterer  on  the  sotirce  axis  has  zero  Doppler  shift,  is  introduced. 
These  three  components  are  added  to  obtain  one  net  Doppler  shift  value  for  each 
cell. 

Subroutine  DOPP  also  serves  a  secondary  purpose.  Included  in  REVMOD  is 
the  capability  to  account  for  Doppler  shift  resulting  from  the  presence  of  current 
layers.  Up  to  eight  layers  are  allowed.  Some  volume  cells  may  fall  on  the  border 
of  two  or  more  current  layers.  In  order  to  obtain  a  better  representation  of  the 
Doppler  for  these  cells,  a  current  split  loop  has  been  set  up  in  DOPP.  This  divides 
a  voltime  cell  into  two  or  more  portions,  depending  on  how  many  current  layers 
that  cell  occupies.  Each  portion  of  the  cell  is  weighted  according  to  the  ratio  of 
its  volume  to  the  total  cell  volume  and  is  processed  individually  for  its  Doppler 
shift  and  frequency  spreading  contributions. 

D.4.3.  CNVOLl 

Once  the  Doppler  shift  has  been  determined  or  a  cell,  the  backscattered 
energy  spectrum  is  shifted  by  the  Doppler  frequency.  This  is  accomplished  by  a 
convolution  of  the  Gatissian  density  distribution  at  the  scatterer  with  an  impulse 
located  at  the  Doppler  shift  value  [27]. 

D.4.4.  SWPDF 

The  backscattered  spectrum  from  the  surface  is  modeled  as  a  combination  of 
two  phenomena.  First  of  all,  a  scattering  function,  due  to  the  roughness  of  the 
surface,  is  computed.  Secondly,  a  surface  wave  contribution  is  introduced  as  a 
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fraction  of  the  total  surface  scattering  contribution.  In  addition  to  the  surface 
scattering  input  variables,  the  surface  wave  speed  and  its  standard  deviation 
are  input.  Thus,  the  surface  reverberation  is  a  combination  of  the  surface 
scattering  component  plus  the  surface  wave  component.  Subroutine  SWPDF 
accounts  for  this  additional  scattering  due  to  surface  waves.  The  probability 
density  function  for  the  surfaM:e  waves  is  added  to  that  for  surface  scatterer. 
The  resulting  composite  probability  density  function  is  normalized  so  that  each 
surface  scattering  cell  is  not  taken  into  account  twice. 

D.4.5.  WTSCAT 

Subroutine  WTSCAT  weights  the  scatterer  speed  distribution  by  the  sonar 
equation  solution  for  each  cell  in  a  horizontal  row.  The  incremental  scattering 
function  is  adjusted  for  a  iHz  band  by  dividing  by  A/,  the  width  of  a  frequency 
bin,  which  is  defined  by  the  sampling  frequency,  /*,  divided  by  the  number  of 
FFT  points,  JV.  The  scattering  functions  for  the  cells  in  an  entire  row  are  then 
acctimulated. 

D.4.6.  RINIT2  /  RINIT3 

Subroutine  RINIT2  is  called  at  the  completion  of  each  rows  of  cells.  The 
vertical  angle,  is  incremented  to  the  center  angle  of  the  next  row  to  be 
processed.  The  azimuthal  angle,  is  reinitialized  to  the  center  angle  of  the 
first  azimuth.  The  surface,  voltime,  and  bottom  components  of  the  scattering 
function  and  total  power  gain  are  separated  for  individual  inspection.  These 
scattering  function  and  power  gain  arrays  are  written  to  a  direct  access  file  for 
plotting  or  as  input  to  other  REVMOD  program  sections.  At  the  end  of  each 
range,  RINIT3  is  called  to  reinitialize  the  arrays.  The  combined  (surface  + 
volume  +  bottom)  scattering  function  and  power  gain  levels  are  written  to  a 
direct  access  file. 
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APPENDIX  E 

INPUT  PARAMETERS 

E.l.  Input  Parameters  for  Figure  3.3 


**•  INPUT  PARAMETERS  »*» 


PARAMETERS  FROM  REVINP  DISK  FILE 
'»»  SOURCE  *** 

SOURCE  SPEED  5.0  M/S  SOURCE  DEPTH  300.0  M 

***  ENVIRONMENT  *** 

SOUND  SPEED  1500.0  M/S  ABSORPTION  COEP.  0.0014  DB/M 
OCEAN  DEPTH  500.0  M 

SURFACE 

SURFACE  SCATTERING  COEFFICIENT  -30.0000  (DB/M«*2) 

STANDARD  DEVIATION  OP  SURF.  SCAT.  SPEED  DIST.  0.1667  M/S 
SURF.  WAVE  DIRECTION  OF  MOTION  REL.  TO  SRC.  HEADING  0.00  DEG 

SURFACE  WAVE  SPEED  0.00  M/S 

STAN.  DEV.  OP  SURF.  WAVE  SPEED  DISTR.  0.167  M/S 

POWER  RATIO  OP  SURF.  WAVE  TO  SURF.  SCAT.  COMP.  1.000 


VOLUME 

VOLUME  SCATTERING  COEFFICIENT  -80.0  DB/M**3 

STANDARD  DEVIATION  OP  VOL.  SCAT.  SPEED  DISTR.  0.17  M/S 


BOTTOM 


BOTTOM  SCATTERING  COEFFICIENT  -35.000  (DB/M**2) 

STANDARD  DEV.  OP  BOTTOM  SCATTERING  SPEED  DISTR.  0.17  M/S 

CURREMTS 


0.0  0.0 

0.0  0.0 

0.0  0.0 


TRANSMITTED  SIGNAL 

CARRIER  FREQUEROr  45000.0  HZ 

PULSE  LENGTH  0.0500000  S 

GRID  CELL  DEPTH  (c*tau/2)  37. SM 


NUMBER  OP  CURRENT 
THICKNESS  (M) 

LAYERS 

0.0 

0 

0.0  0.0 

0.0 

0.0 

0.0 

DIR.  OP  MOT.  (DEG) 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

SPEED  (M/S) 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

BEAMS*** 


TRANSMIT  BEAM  VECTOR  -  THETA  -  0.0  PHI  - 

RECEIVE  BEAM  VECTOR  -  THETA  -  0.0  PHI  - 

BEAM  PATTHIN  FILE  ANGULAR  INCREMENT  3.00  DEG 
NUMBER  BEAM  PATTOIM  POINTS  SPECIFIED  OR  PATTERN  • 


0.0  DEG 
0.0  DEG 
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***  REVERBERATION  MODEL  *** 


BASEBAND  SAMPLING  RATE  BANDWIDTH  1000.0  HZ 

I  OF  BINS  IN  BANDWIDTH  64 

GRID  CELL  WIDTH  IN  AZIMUTH  2.00  DEG 

TOTAL  SPREAD  OP  GRID  AZ.  CEMTHIED  ON  REC.  BEAM  AXIS  10.0  DEG 
(RtID  CELL  HEIGHT  IN  ELEVATION  2.00  DEG 


T-OTAL  SPREAD  OF  GRID  EL.  CENTERED  ON  REC.  BEAM  AXIS  120 


.»»  OUTPUTS  »** 

NUMBER  HP  RANGFS  =  64  RANGE  INCREMENT  =  30. 

RANGE  BEGINNING  =  30.0  M  RANGE  ENDING  = 

SCATTERING  SPECTRUM  VS.  TOTAL  POWER  FLAG  F 
SURFACE, VOLUME, BOTTOM, TOTAL  OR  ALL  FOUR  FLAG  1 


TRANSMIT  BEAM  PATTERN 
O.lOOOE+01  O.lOOOE+01  O.lOOOE+01 
O.lOOOE+01  O.IOOOE+Ol  O.IOOOE+Ol 
O.IOOOE+Ol  O.lOOOE+01  O.IOOOE+Ol 
O.IOOOE+Ol  O.lOOOE+01  O.lOOOE+Ol 
O.lOOOE+01  O.lOOOE+01  Q.IOOOE+Ol 
O.lOOOE+01 


0. lOOOE+01 
O.IOOOE+01 
O.lOOOE+01 
O.lOOOE+01 
O.lOOOE+01 


O.lOOOE+01 

O.lOOOE+01 

O.lOOOE+01 

O.lOOOE+01 

O.lOOOE+01 


RECEIVE  BEAM  PATTERN 
O.lOOOE+01  O.IOOOE+Ol  O.IOOOE+OI 
O.lOOOE+01  O.IOOOE+Ol  O.IOOOE+Ol 
O.lOOOE+01  O.lOOOE+01  O.IOOOE+Ol 
O.lOOOE+01  O.lOOOE+01  O.lOOOE+01 
O.lOOOE+01  O.lOOOE+01  O.lOOOE+01 
O.IOOOE+OI 


O.IOOOE+Ol 

O.lOOOE+01 

O.lOOOE+01 

O.lOOOE+01 

O.lOOOE+01 


O.IOOOE+Ol 

O.lOOOE+01 

O.lOOOE+01 

O.lOOOE+01 

O.lOOOE+01 


.0  DEG 


0  M 

1920.0  M 


O.lOOOE+01 

O.lOOOE+01 

O.lOOOE+01 

O.lOOOE+01 

O.lOOOE+01 


O.lOOOE+01 
0. lOOOE+Ol 
O.lOOOE+01 
O.lOOOE+01 
O.lOOOE+01 


E.2.  Input  Parameters  for  Figure  3.5 


INPUT  PARAMETERS  »»• 


PARAMETERS  FROM  REVINP  DISK  FILE 
SOURCE 

SOURCE  SPEED  5.0  M/S  SOURCE  DEPTH  300.0  M 

•••  ENVIRONMENT 

SOUND  SPEED  1500.0  M/S  ABSORPTION  COEF.  0.0014  DB/M 
OCEAN  DEPTH  500.0  M 

SURFACE 

SURFACE  SCATTERING  COEFFICIENT  -30.0000  (DB/M**2) 

STANDARD  DEVIATION  OF  SURF.  SCAT.  SPEED  DIST.  0.1667  M/S 
SURF.  WAVE  DIRECTION  OF  MOTION  REL.  TO  SRC.  HEADING  0.00  DEG 

SURFACE  WAVE  SPEED  0.00  M/S 

STAN.  DEV.  OP  SURF.  WAVE  SPEED  DISTR.  0.167  M/S 

POWER  RATIO  OP  SURF.  WAVE  TO  SURF.  SCAT.  COMP.  1.000 


VOLUME 

VOLUME  SCATTERING  COEFPICIEOT  -80.0  DB/M»*3 

STANDARD  DEVIATION  OP  VOL.  SCAT.  SPEED  DISTR.  0.17  M/S 

BOTTOM 


BOTTOM  SCATTERING  COEFPICIEMT  -35.000  (DB/M**2) 

STANDARD  DEV.  OF  BOTTOM  SCATTERING  SPEED  DISTR.  0.17  M/S 

CURRENTS 


NUMBER  OF  CURRQiT  LAYERS 
THICKNESS  (M)  0.0 
ILK.  OF  MOT.  (DEG)  0.0 
SPEED  (M/S)  0.0 


0 

0.0  0.0 

0.0  0.0 

0.0  0.0 


0.0  0.0 

0.0  0.0 

0.0  0.0 


0.0  0.0  0.0 

0.0  0.0  0.0 

0.0  0.0  0.0 


TRANSMITTED  SIGNAL 


CARRIER  FRBQUEKCY  45000.0  HZ 

PULSE  LE2IGTH  0.0500000  S 

GRID  CELL  DEPTH  (c*tau/2)  37. SM 


BEAMS*** 


TRANSMIT  BEAM  VECTOR  -  THETA  •  0.0  PHI  - 

RECEIVE  BEAM  VECTOR  -  THETA  -  0.0  PHI  - 

BEAM  PATTERN  FILE  ANGULAR  INCREMBIT  3.00  DEG 
NUMBER  BEAM  PATTERN  POINTS  SPECIFIES  OR  PATTERN  f 


0.0  DEG 
0.0  OEG 
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***  REVERBERATION  MW3EL  *** 


BASmAND  SAMPLING  RATE  BANDWIDTH  1000.0  HZ 

I  OF  BINS  IN  BANDWIDTH  64 

GRID  CELL  WIDTH  IN  AZIMUTH  2.00  DEG 

TOTAL  SPREAD  OF  GRID  AZ.  CENTERED  ON  REC.  BEAM  AXIS 

GRID  CELL  HEIGHT  IN  ELEVATION  3.00  OEG 


10.0  DEG 
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TOTAL  SPREAD  OF  GRID  EL.  CENTERED  ON  REC.  BEAM  AXIS  180. 0  DEG 

OUTPUTS  »»• 

NUMBER  OF  RANGES  =  1  RANGE  INCREMENT  =■  360.0  H 

RANGE  BEGINNING  =  360.0  M  RANGE  E^’DING  =  360.0  M 

SCATTERING  SPECTRUM  VS.  TOTAL  POWER  FLAG  T 
SURFACE, VOLUME, BOTTOM, TOTAL  OR  ALL  FOUR  FLAG 


TRANSMIT  BEAM  PATTERN 
O.lOOOE+01  O.lOOOE+01  O.lOOOE+01 
O.lOOOE+01  O.lOOOE+01  O.lOOOE+01 
O.lOOOE+01  O.lOOOE+01  0 . lOOOE+01 
O.lOOOE+01  O.lOOOE+01  O.lOOOE+01 
O.lOOOE-t-01  O.lOOOE-t-Ol  O.lOOOE+01 
O.lOOOE+01 

RECEIVE  BEAM  PATTERN 
O.lOOOE+01  O.lOOOE+01  O.lOOOE+01 
').1000E+01  O.lOOOE+01  O.IOOOE+Ol 
O.lOOOE+01  O.lOOOE+01  0 . lOOOE+01 
O.IOOOE+Ol  O.lOOOE+01  0 . lOOOE+01 
O.lOOOE+01  O.lOOOE+01  O.lOOOE+01 
O.lOOOE+01 


O.IOOOE+OI  O.lOOOE+01  O.lOOOE+01 
O.lOOOE+01  O.lOOOE+01  O.lOOOE+01 
O.IOOOE+01  O.lOOOE+01  O.lOOOE+01 
O.lOOOE+01  O.lOOOE+01  O.lOOOE+01 
O.lOOOE+01  O.lOOOE+01  O.lOOOE+01 


O.IOOOE+01  O.lOOOE+01  O.lOOOE+01 
O.IOOOE+Ol  O.lOOOE+01  O.lOOOE+01 
O.IOOOE+Ol  O.lOOOE+01  O.lOOOE+01 
O.IOOOE+Ol  O.lOOOE+01  O.lOOOE+01 
O.lOOOE+01  O.lOOOE+01  O.lOOOE+01 


E.3.  Input  Parameters  for  Figure  4.1  and  4.2 


»*•  INPUT  PARAMETERS  *** 


PARAMETERS  FROM  REVINP  DISK  FILE 
»•»  SOURCE  **» 

SOURCE  SPEED  5.0  M/S  SOURCE  DEPTH  300.0  M 

»**  ENVIRONMENT  »*• 

SOUND  SPEED  AT  SOURCE  1500.0  M/S  ABSORPTION  COEF.  VARIES  WITH  SVP 
OCEAN  DEPTH  500.0  M 

NUMBER  OF  BOUNCES  BOUNCES  6  EA.,  SURF.  AND  BOT. 

SURFACE 

SURFACE  SCATTERING  COEFFICIENT  -30.0000  (DB/M**2) 

STANDARD  DEVIATION  OF  SURF.  SCAT.  SPEED  DIST.  0.1667  M/S 
SURF.  WAVE  DIRECTION  OF  MOTION  REL.  TO  VEH.  HEADING  0.00  DEG 

SURFACE  WAVE  SPEED  0.00  M/S 

STAN.  DEV.  OP  SURF.  WAVE  SPEED  DISTR.  0.167  M/S 
POWER  RATIO  OF  SURF.  WAVE  TO  SURF.  SCAT.  COMP.  1.000 
LOSS  PER  SURFACE  BOUNCE  -2.000  dB 


VOLUME 

VOLUME  SCATTERING  COEFFICIENT  -80.0  DB/M**3 

STANDARD  DEVIATION  OP  VOL.  SCAT.  SPEED  DISTR.  0.17  M/S 

BOTTOM 


BOTTOM  SCATTERING  COEFFICIENT  -35.000  (DB/M**25 

STANDARD  DEV.  OF  BOTTOM  SCATTERING  SPEED  DISTR.  0.17  M/S 

LOSS  PER  BOTTOM  BOUNCE  -4.000  dB 


CURRENTS 


NUMBER  OF  CURRERT  LAYERS 
THICKNESS  (M)  0.0 
DIR.  OF  MOT.  (DEG)  0.0 
SPEED  (M/S)  0.0 


0.0 

0.0 

0.0 


0.0 

0.0 

0.0 


0.0 

0.0 

0.0 


0.0 

0.0 

0.0 


0.0  0.0  0.0 

0.0  0.0  0.0 

0.0  0.0  0.0 


TRANSMITTED  SIGNAL 


CARRIER  FREQUENCV  45000.0  HZ 
PULSE  LENGTH  0.. 500000  S 

GRID  CELL  DEPTH  WILL  BE  CALCULATED  IN  REVSCT 


•**  BEAMS*** 


TRANSMIT  BEAM  VECTOR  -  THETA  -  0.0  PHI  - 

RECEIVE  BEAM  VECTOR  -  THETA  -  0.0  PHI  - 

BEAM  PATTERN  FILE  ANGULAR  INCREMENT  3.00  DEG 
NUMBER  BEAM  PATTERN  POINTS  SPECIFIED  OR  PATTERN  I 


0.0  DEG 
0.0  DEG 
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*•*  REVERBERATION  MODEL  *** 


BASEBAND  SAMPLING  RATE  BANDWIDTH  1000.0  HZ 
f  OF  BINS  IN  BANDWIDTH  64 
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GRID  CELL  WIDTH  IN  AZIMUTH  2. 00  DEG 

TOTAL  SPREAD  OF  GRID  AZ.  CENTERED  ON  REC.  BEAM  AXIS  10. ( 
GRID  CELL  HEIGHT  IN  ELEVATION  2.00  DEG 

TOTAL  SPREAD  OF  GRID  EL.  CENTERED  ON  REC.  BEAM  AXIS  180.1 

»**  OUTPUTS  **» 

NUMBER  OF  RANGES  =  64  RANGE  INCREMENT  =  0.040 

RANGE  BEGINNING  =  0.040  S  RANGE  BIDING  =■ 

SCATTERING  SPECTRUM  VS.  TOTAL  POWER  FLAG  F 
SURPACE.VOLUME, BOTTOM, TOTAL  OR  ALL  FOUR  FLAG  1 


TRANSMIT 
O.lOOOE+01 
O.lOOOE+01 
0. lOOOE+01 
0. lOOOE+01 
O.lOOOE+01 
O.lOOOE+01 


BEAM  PATTERN 

O.lOOOE+01 

O.lOOOE+01 

O.lOOOE+01 

O.lOOOE+01 

O.lOOOE+01 


O.lOOOE+01 

O.lOOOE+01 

O.lOOOE+01 

Q.lOOOE+01 

O.IOOOE+Ol 


O.lOOOE+01 

O.lOOOE+01 

O.lOOOE+01 

O.lOOOE+01 

O.lOOOE+01 


0. lOOOE+01 
O.lOOOE+01 
0. lOOOE+01 
O.lOOOE+01 
O.lOOOE+01 


0 

0 

0 

0 

0 


RECEIVE  BEAM  PATTERN 
O.lOOOE+01  O.lOOOE+01  O.lOOOE+01 
O.lOOOE+01  O.lOOOE+01  O.lOOOE+01 
O.lOOOE+01  O.lOOOE-t-01  O.lOOOE+01 
O.lOOOE+01  O.lOOOE+01  O.lOOOE+01 
O.lOOOE+01  O.lOOOE+01  O.lOOOE+01 
O.IOOOEtOI 

NUMDEP  DEPTH  SPEED  ISOCT 


O.lOOOE+01 

O.lOOOE+01 

O.lOOOE+01 

O.lOOOE+01 

O.lOOOE+01 


O.IOOOE+Ol  0 
O.lOOOE+01  0 
O.lOOOE+01  0 
0 . lOOOE+01  0 
O.lOOOE+01  0 


-INV.GRAD.  ALPHA 


1  Sffflrfffl  ISW.'OO  it)6o.M  OB  O.l40fl6-o5 

2  300.00  1500.00  1500.00  0.00  0.1400E-02 

3  0.00  1500.00  0.00  0.00  0.1400B-02 


DEG 

DEG 


S 

1.560  S 


•lOOOE+Ol 
.  lOOOE+01 
-lOOOE+Ol 
■lOOOE+Ol 
.lOOOE+01 


•lOOOE+Ol 

•lOOOE+Ol 

.lOOOE+01 

■lOOOE+Ol 

-lOOOE+Ol 
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E.4.  Input  Parameters  for  Figure  4.3 


INPUT  PARAMETERS  •»» 


PARAMETERS  FROM  REVINP  DISK  FILE 
»»•  SOURCE  •»* 

SOURCE  SPEED  5.0  M/S  SOURCE  DEPTH  100.0  M 

ENVIRONMENT  *»* 

SOUND  SPEED  1500.0  M/S  ABSORPTION  COEF.  0.0014  DB/M 
OCEAN  DEPTH  2000.0  M 

SURFACE 

SURFACE  SCATTERING  COEFFICIENT  -30.0000  (DB/M**2) 

STANDARD  DEVIATION  OF  SURF.  SCAT.  SPEED  DIST.  0.1667  M/S 
SURF.  WAVE  DIRECTION  OF  MOTION  REL.  TO  SRC.  HEADING  O.OD  DEG 

SURFACE  WAVE  SPEED  0.00  M/S 

STAN.  DEV.  OF  SURF.  WAVE  SPEED  DISTR.  0.167  M/S 

POWER  RATIO  OP  SURF.  WAVE  TO  SURF.  SCAT.  COMP.  1.000 


VOLUME 

VOLUME  SCATTERING  COEFFICIENT  -80.0  DB/M**3 

STANDARD  DEVIATION  OF  VOL.  SCAT.  SPEED  DISTR.  0.17  M/S 


BOTTOM 


BOTTOM  SCATTERING  COEFFICIENT  -35.000  (DB/M**2) 

STANDARD  DEV.  OF  BOTTOM  SCATTERING  SPEED  DISTR.  0.17  M/S 

cuRREirrs 


0.0  0.0 

0.0  0.0 

0.0  0.0 


TRANSMITTED  SIGNAL 

CARRIER  FREQUENCY  45000.0  HZ 

PULSE  LENGTH  0.0500000  S 

GRID  CELL  DEPTH  (c*tau/2)  37. SM 


NUMBEN  OF  CURRENT  LAYENS  0 

THICKNESS  (M)  0.0  0.0  0.0  0.0  0.0  0.0 
DIR.  OP  MOT.  (DBGl  0.0  0.0  0.0  0.0  0.0  0.0 
SPEED  (M/S)  0.0  0.0  0.0  0.0  0.0  0.0 


*•*  BEAMS*** 


TRANSMIT  BEAM  VECTOR  -  THETA  -  0.0  PHI  -  0.0  ^ 

RECEIVE  BEAM  VECTOR  -  THETA  -  0.0  PHI  -  0.0  DBG 

BEAM  PATTERN  FILE  ANGULAR  INCNENENT  3.00  DEG 
NUMBER  BEAM  PATTERN  POINTS  SPECIFIED  OR  PATTERN  *  31 


***  REVERBERATION  MODEL  *** 


BASEBAND  SAMPLING  RATE  BANDWIDTH  1000.0  HZ 

I  OF  BINS  IN  BANDWIDTH  64 

GRID  CELL  WIDTH  IN  AZIMUTH  2.00  DEG 

TOTAL  SPREAD  OF  GRID  AZ.  CENTERED  ON  REC.  BEAM  AXIS 

GRID  CELL  HEIGHT  IN  ELEVATION  2.00  DEG 


10.0  DBG 
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TOTAL  SPREAD  OF  GRID  EL.  CENTERED  ON  REC.  3EAM  AXIS  180.0  DEG 


NUMBER  OF  RANGES  =  64  RANGE  INCREMENT  = 

RANGE  BEGINNING  «  96.0  M  RANGE  ENDING 

SCATTERING  SPECTRUM  VS.  TOTAL  POWER  FLAG 
SURFACE, VOLUME, BOTTOM, TOTAL  OR  ALL  FOUR  FLAG 


96.0  M 

6144.0  M 


TRANSMIT  BEAM  PATTERN 
O.lOOOE+01  O.lOOOE+01  O.lOOOE+01 
O.lOOOE+01  O.lOOOE+01  O.lOOOE+01 
O.lOOOE+01  O.lOOOE+01  O.lOOOEi-01 
O.lOOOE+01  O.lOOOE+01  O.lOOOE+01 
O.lOOOE+01  O.lOOOE+01  O.lOOOE+01 
O.lOOOE+01 

RECEIVE  BEAM  PATTERN 
O.lOOOE+01  O.lOOOE+01  O.lOOOE-t-Ol 
O.lOOOE+01  O.lOOOE+01  0 . lOOOE+Ol 
O.lOOOE+01  O.lOOOE+01  O.lOOOE+01 
O.IOOOE+Ol  O.lOOOE+01  O.lOOOE+01 
0.1000E>01  O.lOOOE+01  O.lOOOE+01 
O.lOOOE+01 


O.lOOOE+01  O.lOOOE+01  O.lOOOE+01 
O.IOOOE+Ol  O.lOOOE+01  O.lOOOE+01 
O.lOOOE+01  O.lOOOE+01  O.lOOOE+01 
O.IOOOE+Ol  O.lOOOE+01  O.lOOOE+01 
O.lOOOE+01  O.lOOOE+01  O.lOOOE+01 


O.lOOOE+01  O.lOOOE+01  O.lOOOE+01 
O.lOOOE+01  O.lOOOE+01  O.lOOOE+01 
O.lOOOE+01  O.lOOOE+01  O.IOOOE+Ol 
O.lOOOE+01  O.lOOOE+01  C.lOOOE+01 
O.IOOOE+Ol  O.lOOOE+01  O.lOOOE+01 


E.5.  Input  Parameters  for  Figure  4.5 


»*»  INPUT  PARAMETERS  *** 


PARAMETERS  FROM  REVINP  DISK  FILE 
»»•  SOURCE  »•» 

SOURCE  SPEED  5.0  M/S  SOURCE  DEPTH  100.0  M 

»*•  ENVIRONMENT  **• 

SOUND  SPEED  AT  SOURCE  1510.0  M/S  ABSORPTION  COEF.  VARIES  WITH  SVp 
I'CEAN  DEPTH  <:000.0  M 

NUMBER  OF  BOUNCES  BOUNCES  6  EA..  SURF.  AND  BOT. 

SURFACE 

SURFACE  SCATTERING  COEFFICIENT  -30.0000  (DB/M**2) 

STANDARD  DEVIATION  OF  SURF.  SCAT.  SPEED  DIST.  0.1667  M/S 
SURF.  WAVE  DIRECTION  OF  MOTION  REL.  TO  VEH.  HEADING  0.00  DEG 

SURFACE  WAVE  SPEED  0.00  M/S 

STAN.  DEV.  OP  SURF.  WAVE  SPEED  DISTR.  0.167  M/S 
POWER  RATIO  OP  SURF.  WAVE  TO  SURF.  SCAT.  COMP.  1.000 
LOSS  PER  SURFACE  BOUNCE  -2.000  dB 


VOLUME 

VOLUME  SCATTERING  COEFFICIENT  -80.0  DB/M**3 

STANDARD  DEVIATION  OF  VOL.  SCAT.  SPEED  DISTR.  0.17  M/S 


BOTTOM 


BOTTOM  SCATTERING  COEFFICIENT  -35.000  {DB/M**2) 

STANDARD  DEV.  OF  BOTTOM  SCATTERING  SPEED  DISTR.  0.17  M/S 

LOSS  PER  BOTTOM  BOUNCE  -4.000  dB 

CURRENTS 


NUMBER  OP  CURRENT  LAYERS 
THICKNESS  (M)  0.0 
DIR.  OP  MOT.  (DBG)  0.0 
SPEED  (M/S)  0.0 


0.0 

0.0 

0.0 


0.0 

0.0 

0.0 


0.0 

0.0 

0.0 


0.0 

0.0 

0.0 


0.0 

0.0 

0.0 


0.0 

0.0 

0.0 


0.0 

0.0 

0.0 


TRANSMITTED  SIGNAL 


CARRIER  FREQUENCY  45000.0  HZ 
PULSE  LENGTH  0.0500000  S 

GRID  CELL  DEPTH  WILL  BE  CALCULATED  IN  REVSCTT 


*•*  BEAMS*** 

TRANSMIT  BEAM  VECTOR  -  THETTA  -  0.0  PHI  -  0.0  DEG 

RECEIVE  BEAM  VECTOR  -  THETA  -  0.0  PHI  -  0.0  DEG 

BEAM  PATTERN  PILE  ANGULAR  INCREMENT  3.00  DEG 
NUMBS)  BEAM  PATTERN  POINTS  SPECIFIED  OR  PATTERN  *  31 


***  REVERBERATION  MODEL  *** 

BASEBAND  SAMPLING  RATE  BANDWIDTH  1000.0  HZ 
t  OF  BINS  IN  BANDWIDTH  64 
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GRID  CELL  WIDTH  IH  AZIMUTH  2.00  DEG 

TOTAL  SPREAD  OF  GRID  AZ.  CENTERED  ON  REC.  BEAM  AXIS  10. 

GRID  CELL  HEIGHT  IN  ELEVATION  2.00  DEG 

total  spread  of  grid  el.  centered  on  REC.  BEAM  AXIS  180. 

»»»  OUTPUTS  *»» 

NUMBER  OF  RANGES  =  64  RANGE  INCREMENT  »  0.128 

RANGE  BEGINNING  =  0.128  S  RANGE  ENDING  = 

SCATTERING  SPECTRUM  VS.  TOTAL  POVfER  FLAG  T 
SURFACE, VOLUME, BOTTOM, TOTAL  OR  ALL  FOUR  FLAG  1 


TRANSMIT 

BEAM  PATTERN 

0. lOOOE+Ol 

0. lOOOE+01 

O.IOOOE+Ol 

O.lOOOE+01 

O.lOOOE+01 

0 

O.IOOOE+Ol 

0. lOOOE+01 

O.lOOOE+01 

O.lOOOE+01 

O.lOOOE+01 

0 

0. lOOOE+01 

0. lOOOE+01 

O.lOOOE+01 

O.lOOOE+01 

O.lOOOE+01 

0 

0. lOOOE+01 

0. lOOOE+01 

O.lOOOE+01 

O.lOOOE+01 

O.lOOOE+01 

0 

0 . lOOOE+Ol 
O.lOOOE+Ol 

0, lOOOE+01 

O.IOOOE+Ol 

O.lOOOE+01 

O.lOOOE+01 

0 

RECEIVE 

BEAM  PATTERN 

O.lQOOE+01 

O.lOOOE+01 

O.IOOOE+Ol 

O.lOOOE+01 
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APPENDIX  F 

CONSTANT  FACTOR  IN  DOPPLER  EQUATIONS 
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The  Doppler  shift  equation  for  a  moving  source  has  been  defined  as 

/£>  =  /o  ^  cos  7^  ,  (F.l) 

where  /<>  is  the  original  frequency,  v  is  the  source  velocity,  c  is  the  speed  of  sound, 
and  7  is  the  angle  at  which  the  Doppler  is  realized.  In  an  isovelocity  environment, 
the  incident  angle,  7^,  at  a  scatterer,  is  identical  to  the  angle,  7,  at  the  source. 
Therefore,  the  Doppler  equation  for  a  moving  scatterer  is 

/d  =  /o  ^  cos  7^  .  (F,2) 

The  difference  in  the  Doppler-shifted  frequency  from  the  original  is 


for  source  motion,  and 


A/  = - cos  7 


A  /  2/ott 

A/  * - cos  7 


for  scatterer  motion.  To  simplify  the  calculations,  the  constant  factor, 

2/0 

-—cos  7, 


(F.3) 


(FA) 


(F.5) 


is  used  in  program  SCATTER  for  Doppler  due  to  either  source  or  scatterer 
motion. 


When  a  non-isovelocity  environment  is  present,  this  simplification  might  not 
seem  valid,  although  it  is.  In  general,  in  a  non-isovelocity  scenario,  the  angle,  7, 
at  the  source  is  not  equal  to  the  auigle,  7s,  at  the  scatterer.  Nor  is  the  speed  of 
sound  at  the  source  equal  to  the  speed  of  sound  at  the  scatterer.  However,  there 
is  a  relation  between  the  angles  and  speeds  of  sound  from  one  point  to  another. 
This  relation  has  been  introduced  as  Snell’s  law,  and  it  states  that 

siny  s\n'^g 
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where  c  is  the  speed  of  sound  at  the  source  and  is  the  speed  of  sound  at  the 
scatterer.  The  angles,  7'  and  7^,  are  the  angles  at  the  source  and  the  scatterer, 
respectively.  These  angles,  however,  must  be  measured  with  the  convention  of 
zero  degrees  being  straight  downward,  so  that  Snell’s  law  applies.  That  is  to  say 
that  zero  degrees  must  be  perpendicular  to  the  horizontally  stratified  layers.  To 
convert  the  angles,  'y'  and  7',  to  the  actual  angles  in  SCATTER,  7  and  7,,  a 
shift  of  ninety  degrees  must  be  introduced.  Snell’s  law  may  now  be  expressed  as 


sin7  +  90  _  sin 7a  +  90 

C  Ca 


(F.l) 


This  simplifies  to 


cos  7  cos  7a 


(E.8) 


C  Ca 

Therefore,  when  the  Doppler  shift  at  a  scattering  element  is  determined,  the 
factor 

2/0 


cos  7, 


(F.9) 


is  a  constant,  whether  it  is  measured  at  the  soTirce  or  at  the  scatterer. 


APPENDIX  H 

RAYTRACING  SUBROUTINE  FLOWCHARTS 


H.l.  CRITICAL  ANGLES 


143 


H.2.  RAYSET 


C  ENTER  ) 


INPUT:  DATA  FROM 
SUBROUTINE  CRITICAL- 
ANGLES,  RANGE,  PULSE 
LENGTH 


GRID  CELL  DEPTH  (GC0)=l/2  PULSE  LENGTH 

FRONT  OF  SCATTERING 
VOLUMEsRANGE 


BACK  OF  SCATTERING 
VOLUME=RANGE+GCD 

DETERMINE,  FROM  THE 
CRITICAL  ANGLES,  HOW 
MANY  SECTORS  TO  BE  SEARCHED 


TIME:  SURFACE 
2*^0  TijyflE;  bottom 


(  RETURN  ) 
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H.5.  BROKEN 


ENTER 


DETERMINE,  THROUGH  THE  USE 
OF  dt/d(b,  WHICH  ANGLE 
BOUNDING  THE  SECTOR  ENSONIFIES 
THE  BOUNDARY  WITHIN  THE 
SCATTERING  SHELL.  (THIS  IS 
ONE  SOLUTION  ANGLE.) 


ITERATE  FOR  THE  OTHER 
SOLUTION  ANGLE 


ITERATION 

SECTION 


RETURN 


h.6.  iteration  section 


CUT  SECTOR  IN  HALF 
TO  OBTAIN  TRUL 
SOLUTION  ANGLE 


<(^TIMEBOUND^ 


COMBINE  PARTIAL  TRAVEL  TIMES 
FROM  TIMEBOUND  WITH  STATE 
NUMBERS  K.  KA.  AND  KB 
TO  OBTAIN  A  TRML 
SOLUTION  TIME 


COMPARE  TRIAL  SOLUTION 
TIME  WITH  THE  TIME  AT 
ONE  END  OF  THE  SHELL 


USE  TO 
DETERMINE  WHICH  HALF 
OF  THE  PREVIOUS 
SECTOR  TO  SEARCH 
NEXT 


CUT  THIS  NEW  SECTOR 
•N  half  to  OBTAIN  A 
NEW  TRIAL  ANGLE 


OIMEBOUNoS 

COMPARE  NEWi 
TRUL  TIME 
WITH  OLD 
TRUL  TIME 


I  CALCULATE  THE  LENGTH 
I  OP  THE  ENSONinED  PATCH 

I 


(  RETURN  ^ 
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H.7.  SOLBND 
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H.8.  VOLANGS 


ORDER  THIS  COMPLETE  SET  OF 
'*VOLANGS"  FROM  0  TO  180  DEGREES 


H.9.  SURSEQ 


H.ll.  VOLSEQ 


I  OITtllMINC  THC  LASGfST  GSIO  CCU  HftONT  I 

)  WHICH  IS  less  THAN  Off  iOUAi  TO  THi  INffUT  I 

I  MAXIMUM  CEIL  HSIGHT.  SUCH  THAT  THE  SOLUTION  I 


53 


H.12.  BOTSEQ 


(^EMTEW  ) 
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H.13.  RAYDOP 


GizfD 
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H.15.  TURN 


H.16.  pieces 
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C  ENTER 


_ I- 

CONVERT  ALL 
ANGLES  FROM 
DEGREES  TO 
RADIANS 


_ 1 _ 

DETERMINE  THE 
ARCLENGTH  THROUGH 
THE  PRESENT  LAYER 


DETERMINE  THE 
ABSORPTION  LOSS 
FOR  THIS  LAYER 


DETERMINE  THE 
HORIZONTAL  DISTANCE 
TRAVELED  THROUGH 
THIS  LAYER 


RETURN 
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